Python: taking the GLCM of a non-rectangular region

Question

I have been using the SLIC implementation of skimage to segment images in superpixels. I would like to use GLCMs to extract additional features from these superpixels for a classification problem. These superpixels are not rectangular. In MATLAB you can set pixels to NaN and they will be ignored by the algorithm (link). I could use this to make bounding boxes around the superpixels and then just setting the unused pixels to NaN.

The greycomatrix function in skimage does not work entirely the same as the MATLAB implementation however. When setting pixels to NaN the function fails on an assert to check if all values are larger than 0.

Is there a Python implementation available which is able to work with nonrectangular ROIs?

Answer 1

Although mahotas is also an excellent computer vision library, there's no need to stop using skimage to do this.

What is necessary, as @Tonechas has pointed out, is to set those NaN values to an integer, since np.nan has type float and the greycomatrix function requires an array of integers.

The easiest option would be setting those NaN's to zero but, if you already have zero values in your pixels and don't want to mix them, you can choose any other constant. After that, all you have to do is filter that chosen value (once again, generally zero) out of the GLCM.

To understand what this means, let's see what skimage tells us about the output of the greycomatrix function:

4-D ndarray

[...] The value P[i,j,d,theta] is the number of times that grey-level j occurs at a distance d and at an angle theta from grey-level i. If normed is False, the output is of type uint32, otherwise it is float64. The dimensions are: levels x levels x number of distances x number of angles.

In other words, the first two dimensions of the array define a matrix that tells us how many times two different values are a certain distant apart. Note that the GLCM does not keep the shape of the input array. Those rows and columns are telling us how the values relate.

Knowing this, it's easy to filter out the values outside our ROI (imagine we've set those NaN's to zero):

glcm = greycomatrix(img, [1], [0])  # Calculate the GLCM "one pixel to the right"
filt_glcm = glcm[1:, 1:, :, :]           # Filter out the first row and column

Now you could easily calculate the Haralick properties of your filtered GLCM. For example:

greycoprops(filt_glcm, prop='contrast')

Answer 2

The issue is you have to pass an integer array to greycomatrix, but np.nan has type float (take a look at this thread for details). As a result you cannot encode the pixels outside the ROI as NaN.

An approximate workaround for dealing with non-rectangular ROI's would be setting the pixels outside the ROI to 0 and using the function haralick from mahotas library. This function returns 13 Haralick features extracted from four different GLCM's, corresponding to the four 2-D orientations and the particular value of the distance parameter.

From the documentation:

ignore_zeros can be used to have the function ignore any zero-valued pixels (as background).

In summary, you need to mask those pixels that fall outside the ROI, and set ignore_zeros to True in the call to haralick.

---

DEMO

To begin with, let us generate some mock data:

In [213]: import numpy as np

In [214]: shape = (3, 4)

In [215]: levels = 8

In [216]: np.random.seed(2017)

In [217]: x = np.random.randint(0, levels, size=shape)

In [218]: x
Out[218]: 
array([[3, 1, 6, 5],
       [2, 0, 2, 2],
       [3, 7, 7, 7]])

Then we have to remove all the zeros from the image as in this approach the zero intensity level is reserved for the pixels outside the ROI. It is worth to point out that merging intensities 0 and 1 into a single intensity 1 introduces inaccuracies in the results.

In [219]: x[x == 0] = 1

In [220]: x
Out[220]: 
array([[3, 1, 6, 5],
       [2, 1, 2, 2],
       [3, 7, 7, 7]])

Next step consists in defining a mask for the pixels outside the ROI (in this toy example, the four corners of the image) and setting those pixels to 0.

In [221]: non_roi = np.zeros(shape=shape, dtype=np.bool)

In [222]: non_roi[np.ix_([0, -1], [0, -1])] = True

In [223]: non_roi
Out[223]: 
array([[ True, False, False,  True],
       [False, False, False, False],
       [ True, False, False,  True]], dtype=bool)

In [224]: x[non_roi] = 0

In [225]: x
Out[225]: 
array([[0, 1, 6, 0],
       [2, 1, 2, 2],
       [0, 7, 7, 0]])

We can now perform feature extraction from the GLCM's of a non-rectangular ROI:

In [226]: import mahotas.features.texture as mht

In [227]: features = mht.haralick(x, ignore_zeros=True)

In [228]: features.size
Out[228]: 52

In [229]: features.ravel()
Out[229]: array([ 0.18      ,  5.4       ,  0.5254833 , ...,  0.81127812,
       -0.68810414,  0.96300727])

It may be useful to check how the co-occurrence matrices look. For example, the "pixel to the right" GLCM would be:

In [230]: mht.cooccurence(x, 0)
Out[230]: 
array([[0, 1, 0, ..., 0, 1, 2],
       [1, 0, 2, ..., 0, 1, 0],
       [0, 2, 2, ..., 0, 0, 0],
       ..., 
       [0, 0, 0, ..., 0, 0, 0],
       [1, 1, 0, ..., 0, 0, 0],
       [2, 0, 0, ..., 0, 0, 2]])