R: Calculate Grey-Level-Co-Occurence-Matrix (GLCM) for an image

Question

To investigate on the distribution of pixelvalues in an image, I want to compute a Grey-Level-Co-Occurence-Matrix (GLCM) for entire Images (NO sliding/moving Windows). The idea is to receive a single value (for "mean", "variance", "homogeneity", "contrast", "dissimilarity", "entropy", "second_moment", "correlation") for every image, to compare the images among each other regarding their distribution of pixelvalues.

e.g.:

image 1:

0 0 0 0 
0 0 0 1
0 0 1 1
0 1 1 1

image 2:

1 0 0 1 
0 1 0 0
0 0 1 0
1 0 0 1

image 3:

1 1 1 0
1 1 0 0
1 0 0 0
0 0 0 0

All of These 3 images have got the same statistics (mean, max, min, …), nevertheless the distribution of the pixelvalues is completely different. To find kind of a measure to describe that difference, I want to compute the GLCM´s for each of these images.

I am using the package "glcm" so far, a fantastic package for texture-analysis by Alex Zvoleff. Unfortunately it´s just possible to use it with a sliding/moving window… But since I want to receive one single value for every image per statistical measure it seems to be useless for me... Is there anyone who can help an R-Rookie like me out with that? :)

install.packages("glcm")
library(glcm)
# install and load package "glcm"
# see URL:http://azvoleff.com/articles/calculating-image-textures-with-glcm/

values <- seq(1, c(12*12), 1)
values_mtx <- matrix(data = values, nrow = 12, ncol = 12, byrow = TRUE)
# create an "image"

values_mtx_small <- values_mtx[-12, -12]
# since you have to use a sliding/moving window in glcm::glcm() give the image # ...an odd number of rows and cols by deleting the last row and last column

values_raster_small <- raster(values_mtx_small)
# create rasterlayer-object

values_textures <- glcm::glcm(values_raster_small, window = c((nrow(values_raster_small)-2), (ncol(values_raster_small)-2)), shift=list(c(0,1), c(1,1), c(1,0), c(1,-1)), statistics = c("mean", "variance", "homogeneity", "contrast", "dissimilarity", "entropy", "second_moment", "correlation"), min_x = NULL, max_x = NULL, na_opt = "ignore", na_val = NA, asinteger = FALSE)
# compute a GLCM for the image with a maximum size for the moving window to
# ...receive a "measure" for the image

values_textures_mean <- as.matrix(values_textures$glcm_mean)
# extract the calculated GLCM_mean data

values_textures_mean
# get an Output

   [,1] [,2] [,3] [,4] [,5]      [,6] [,7] [,8] [,9] [,10] [,11]
 [1,]   NA   NA   NA   NA   NA        NA   NA   NA   NA    NA    NA
 [2,]   NA   NA   NA   NA   NA        NA   NA   NA   NA    NA    NA
 [3,]   NA   NA   NA   NA   NA        NA   NA   NA   NA    NA    NA
 [4,]   NA   NA   NA   NA   NA        NA   NA   NA   NA    NA    NA
 [5,]   NA   NA   NA   NA   NA 0.4589603   NA   NA   NA    NA    NA
 [6,]   NA   NA   NA   NA   NA 0.5516493   NA   NA   NA    NA    NA
 [7,]   NA   NA   NA   NA   NA        NA   NA   NA   NA    NA    NA
 [8,]   NA   NA   NA   NA   NA        NA   NA   NA   NA    NA    NA
 [9,]   NA   NA   NA   NA   NA        NA   NA   NA   NA    NA    NA
[10,]   NA   NA   NA   NA   NA        NA   NA   NA   NA    NA    NA
[11,]   NA   NA   NA   NA   NA        NA   NA   NA   NA    NA    NA
# unfortunately two numbers as "measure" are left…

Answer 1

This suggestion might provided the tools needed to get at the answer through the package EBImage. The complete answer would likely require applying additional data reduction techniques and statistical analysis to the results from the textural analysis demonstrated here.

# EBImage needed through Bioconductor, which uses BiocManager
  if (!require(EBImage)) {
    if (!requireNamespace("BiocManager", quietly = TRUE))
      install.packages("BiocManager")
    BiocManager::install("EBImage")
    library(EBImage)
  }

For EBImage, a binary mask is required to define objects for subsequent analysis. In this case, the entire image (array) seems to serve as the object of analysis so a binary mask covering the entire image is created and then modified to replicate the example.

# Create three 32 x 32 images similar to the example
  mask <- Image(1, dim = c(32, 32))
  img1 <- img2 <- img3 <- mask
  img1[upper.tri(img1)] <- 0
  nzero <- sum(img1 == 0)
  img2[sample(32*32, nzero)] <- 0
  img3[lower.tri(img3)] <- 0

# Combine three images into a single 64 x 64 x 3 array for simplicity
  img <- combine(img1, img2, img3)

# Verify similarity of global properties of each image
  apply(img, 3, mean)
> [1] 0.515625 0.515625 0.515625
  apply(img, 3, sd)
> [1] 0.5 0.5 0.5

Haralick features computes rotational invariant textural properties from the gray-level co-occurrence matrix. The parameter haralick.scales is used to specify the expected repeating scale for the textural patterns. The default uses c(1, 2) to look for repeats every 1 and 2 pixels. Here I just limit it to 1 pixel.

I have to admit that I use it without fully understanding it. One helpful resource may be a post by Earl Glynn. Also, a question answered on the Bioconductor about computing Haralick features provides great information that's hard to find.

# Introduce and apply the computeFeatures.haralick function at a scale of 1
# The first line simply captures the names and properties of the features
  props <- computeFeatures.haralick(properties = TRUE, haralick.scales = 1)

# Apply computeFeatures.haralick to each of the 3 dimensions (frames)
  m <- sapply(getFrames(img),
    function(ref) computeFeatures.haralick(mask, ref, haralick.scales = 1))

# Add meaningful row and column names to the resulting matrix
  rownames(m) <- props$name
  colnames(m) <- paste0("img", 1:3)
  print(round(m, 4))
>               img1      img2      img3
> h.asm.s1    0.4702    0.2504    0.4692
> h.con.s1   30.7417  480.7583   30.7417
> h.cor.s1    0.9359   -0.0013    0.9360
> h.var.s1  240.6937  241.0729  241.1896
> h.idm.s1    0.9680    0.5003    0.9680
> h.sav.s1   34.4917   33.8417   33.4917
> h.sva.s1 2093.5247 1594.4603 2028.1987
> h.sen.s1    0.3524    0.4511    0.3528
> h.ent.s1    0.3620    0.6017    0.3625
> h.dva.s1    0.0000    0.0000    0.0000
> h.den.s1    0.0137    0.1506    0.0137
> h.f12.s1    0.7954    0.0000    0.7957
> h.f13.s1    0.6165    0.0008    0.6169

Here I use a heatmap to visualize and organize the 13 Haralick parameters. The plot pretty clearly shows that images 1 and 3 are rather similar and quite different from image 2. Still, differences between image 1 and 3 can be seen.

The matrix used for this heatmap, especially if it was generated from many more images, could be scaled and further analyzed by principle components analysis to identify related images.

  heatmap(m)

To learn more about EBImage see the the package vignette.

Answer 2

My R package GLCMTextures is mainly meant to deal with spatial raster data like glcm, but it should be able to do this too. You'll have to tabulate a GLCM for each of the four shifts [c(1, 0), c(1, 1), c(0, 1), c(-1, 1)] individually and then average the texture metrics of each type to get directionally invariant measures.

library(GLCMTextures)
library(raster)

# create an "image"
values_mtx <- matrix(data = seq(1, c(12*12), 1), nrow = 12, ncol = 12, byrow = TRUE)
values_mtx_raster<- raster(values_mtx) #Make it a raster
values_mtx_raster_quantized<- quantize_raster(values_mtx_raster, n_levels = 32, method = "equal prob") #make values integers from 0-31

plot(values_mtx_raster_quantized)
text(values_mtx_raster_quantized)

values_mtx_quantized<- as.matrix(values_mtx_raster_quantized) #Make it a matrix

glcm_10<- make_glcm(values_mtx_quantized, n_levels = 32, shift = c(1,0), na.rm = FALSE, normalize = TRUE) #tabulate glcm with xshift=1, yshift=0 (i.e. pixel to the right)
glcm_metrics(glcm_10)
# glcm_contrast glcm_dissimilarity   glcm_homogeneity           glcm_ASM       glcm_entropy          glcm_mean      glcm_variance   glcm_correlation 
# 0.21212121         0.21212121         0.89393939         0.02100551         4.08422180        15.50000000        84.25000000         0.99874112