How to calculate the grey-level co-occurrence matrix with only some pixels?

Question

I want to calculate the skimage.feature.texture.greycomatrix(). Now, I have tried it with the following parameters:

skimage.feature.texture.greycomatrix(image, [1], [0], levels=256, symmetric=False, normed=False)

But, what I want to do, is to calculate this matrix with only some pixels.

I have this mask:

And, I have this image:

And, I only want to calculate the intersection pixels between the white mask and the picture. I don't want to calculate the mask with black pixels or other elements left over, only with the pixels of the image that intersect with the white area of the mask like this:

How can I do this?

Answer 1

Solely using skimage.measure.greycomatrix, I don't think you'll be able to achieve that, since there's no mask parameter or similar. So, I guess, you'd need to build something on your own here.

My idea would be to set all unwanted pixels to -1, get the desired correspondences with respect to the (x, y) offset, filter out all those containing at least one -1 entry, and count the occurences of all remaining correspondences. Such, you neglect all correspondences which contain a pixel outside the mask.

To show that for a simple example (cf. the first example from skimage.measure.greycomatrix):

import numpy as np

# https://scikit-image.org/docs/dev/api/skimage.feature.html#greycomatrix
image = np.array([[0, 0, 1, 1],
                  [0, 0, 1, 1],
                  [0, 2, 2, 2],
                  [2, 2, 3, 3]], dtype=np.uint8)
h, w = image.shape

# Set offset and levels (assuming range [0 ... max] here)
x_off, y_off = (1, 0)
levels = np.max(image) + 1

# Temporary matrix
temp = -np.ones((h+y_off, w+x_off), int)
temp[:h, :w] = image
temp = temp[y_off:, x_off:]

# Get correspondences
temp = np.vstack([image.flatten(), temp.flatten()]).T

# Filter correspondences with at least one entry -1, and count correspondences
temp = temp[~np.any(temp == -1, axis=1), :]
temp = np.unique(temp, axis=0, return_counts=True)

# Create GLCM
glcm = np.zeros((levels, levels), int)
glcm[temp[0][:, 0], temp[0][:, 1]] = temp[1]
print(glcm)
# [[2 2 1 0]
#  [0 2 0 0]
#  [0 0 3 1]
#  [0 0 0 1]]

That's in line with the results of skimage.measure.greycomatrix.

For your actual image and mask, that'd be the full code (sorry for using OpenCV here; for me, it's somehow faster to prototype):

import cv2
import matplotlib.pyplot as plt
import numpy as np

# Read image and mask
image = cv2.imread('us.png', cv2.IMREAD_GRAYSCALE)[:857, :1211]
mask = cv2.imread('mask.png', cv2.IMREAD_GRAYSCALE)[:857, :1211]

# Get bounding rectangle, crop image and mask
x, y, w, h = cv2.boundingRect(mask)
image = image[y:y+h, x:x+w]
mask = mask[y:y+h, x:x+w]

# DEBUG: Count non-zero pixels of mask
print('# non-zero pixels mask:', cv2.countNonZero(mask))

# Set all non-masked pixels in image to -1
image = image.astype(int)
image[~mask.astype(bool)] = -1
h, w = image.shape

# Set offset and levels (assuming range [0 ... max] here)
x_off, y_off = (1, 0)
levels = np.max(image) + 1

# Temporary matrix
temp = -np.ones((h+y_off, w+x_off), int)
temp[:h, :w] = image
temp = temp[y_off:, x_off:]

# Get correspondences
temp = np.vstack([image.flatten(), temp.flatten()]).T

# Filter correspondences with at least one entry -1, and count correspondences
temp = temp[~np.any(temp == -1, axis=1), :]
temp = np.unique(temp, axis=0, return_counts=True)

# Create GLCM
glcm = np.zeros((levels, levels), int)
glcm[temp[0][:, 0], temp[0][:, 1]] = temp[1]

# DEBUG: Count number of correspondences found in GLCM
print('# correspondences found in GLCM:', np.sum(glcm))

# Some visualization
plt.figure(1, figsize=(8, 8))
plt.imshow(glcm, cmap='hot'), plt.colorbar()
plt.tight_layout(), plt.show()

The GLCM would be too large to print, so the plot:

Let's see the additional DEBUG output:

# non-zero pixels mask: 18636
# correspondences found in GLCM: 18259

There are 18,636 non-zero pixels in your mask, but only 18,259 correspondences were counted in the GLCM. I haven't counted those pixels to be neglected, but for offset (1, 0), the difference of about 400 seems reasonable.

----------------------------------------
System information
----------------------------------------
Platform:      Windows-10-10.0.16299-SP0
Python:        3.9.1
PyCharm:       2021.1.1
Matplotlib:    3.4.1
NumPy:         1.20.2
OpenCV:        4.5.1
----------------------------------------