2

I'm trying to calculate some textural measures using the GLCM described by Haralick (energy, homogeneity, etc.) for a series of 4 band (R, G, B, NIR) aerial photographs that I have. I have tried this on a subset but am ending up with an image that is mostly blank. My current understanding is that it has to do with the greyscaling and the levels parameter but I can't figure it out.

My date is very large (several GB) so I'm trying to be efficient by using the module RIOS (reads an image in as a 400×400×nbands numpy array, processes the data and writes out to an output image).

My input scene can be found here (200 MB).

My output image looks like (this may be difficult to see as the black pixels are very small):

output

My code is:

#Set up input and output filenames
infiles = applier.FilenameAssociations()
infiles.image1 = "infile.tif"

outfiles = applier.FilenameAssociations()
outfiles.outimage = "outfile.tif"

controls = applier.ApplierControls()
controls.progress = cuiprogress.CUIProgressBar()
# I ultimately want to use a window here
# which RIOS easily allows you to set up.
# For a 3x3 the overlap is 1, 5x5 overlap is 2 etc
#controls.setOverlap(4)

def doFilter(info, infiles, outfiles, controls=controls):
    grayImg = img_as_ubyte(color.rgb2gray(infiles.image1[3]))
    g = greycomatrix(grayImg, distances=[1], angles=[0, np.pi/4, np.pi/2, 3*np.pi/4], symmetric=True, normed=True)
    filtered = greycoprops(g, 'energy')
    # create 3d image from 2d array
    outfiles.outimage = numpy.expand_dims(filtered, axis=0)


applier.apply(doFilter, infiles, outfiles, controls=controls)

Obviously there is something wrong here as my output is not as I expect. My guess that it is to do with the 'levels' parameter. I have been pointed to an explanation here: Black line in GLCM result which explains the parameter well but I am unable to improve my result.

Can someone explain to me why my result is coming out as shown and how I can remedy it?

Community
  • 1
  • 1
GeoMonkey
  • 1,615
  • 7
  • 28
  • 56
  • Your image is binary, all the pixel intensities are either `0`or `255`. Execute `np.unique()` to convince yourself. The GLCM for such an image would have only four non-zero entries. – Tonechas Apr 06 '17 at 22:30
  • numpy.unique yields [ 21 22 23 24 25 26 27 28 29 30 3......... 186 187 188 189 190 191 192 193 194 195 196 197 198] – GeoMonkey Apr 06 '17 at 22:39
  • I ran this code: `import numpy as np` `from skimage import io` `x = io.imread('https://i.stack.imgur.com/EyCI1.png')` `np.unique(x)` and obtained: `array([ 0, 255], dtype=uint8)` – Tonechas Apr 06 '17 at 22:45
  • ah, my bad, thats a screenshot of my output. I will add my real data in – GeoMonkey Apr 06 '17 at 22:50

1 Answers1

2

The code below computes the GLCM corresponding to an offset "1-pixel offset upwards" from the NIR band of your tif image:

import numpy as np
from skimage import io
from skimage.feature import greycomatrix, greycoprops

x = io.imread('m_2909112_se_15_1_20150826.tif')
nir = x[:, :, 3]

glcm = greycomatrix(nir, [1], [np.pi/2], levels=256, normed=True, symmetric=True)

This is how nir looks:

NIR band

The effect of setting the parameter normed to True is that the computed GLCM is divided by its total sum, and as a result the elements of glcm have rather small values. Here's a sample:

In [48]: np.set_printoptions(precision=3)

In [49]: glcm[:5, :5, 0, 0]
Out[49]: 
array([[  0.000e+00,   0.000e+00,   0.000e+00,   0.000e+00,   0.000e+00],
       [  0.000e+00,   2.725e-03,   6.940e-05,   3.725e-05,   2.426e-05],
       [  0.000e+00,   6.940e-05,   1.709e-04,   4.103e-05,   2.216e-05],
       [  0.000e+00,   3.725e-05,   4.103e-05,   4.311e-04,   4.222e-05],
       [  0.000e+00,   2.426e-05,   2.216e-05,   4.222e-05,   5.972e-05]])

To display glcm as an image you need to rescale it, for example like this:

from skimage.exposure import rescale_intensity
scaled = rescale_intensity(glcm[:,:,0,0])
io.imshow(scaled)

scaled

Tonechas
  • 13,398
  • 16
  • 46
  • 80
  • I have implemented what you suggest but my result is the same. Ive looked into the outputs a bit more using print statements. the data is being read in in blocks as I expect, like: [174 169 176 179 200 203 200 183 175 192 186 180 183 181 187 178 170 177 171 ......173 166 155 148 166 162 138 115 146 163 172 170 172 173 173 180 181 173 171 171] but the results of the greycomatrix are often empty, like [0 0 0 0] with only a few that are populated, like [0 1 0 1]. This is essentially passing blank arrays to the greycoprops command. – GeoMonkey Apr 06 '17 at 23:54