R:计算图像的灰度共生矩阵(GLCM)

Question

为了研究图像中像素值的分布，我想计算整个图像(没有滑动/移动窗口)的灰度共生矩阵(GLCM)。其思想是为每个图像接收单个值(代表“均值”、“方差”、“同质性”、“对比度”、“相异度”、“熵”、"second_moment“、”相关性“)，以比较图像之间的像素值分布。

例如：

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

所有这3个图像都有相同的统计数据(均值，最大值，最小值，…)然而，像素值的分布是完全不同的。为了找到一种度量来描述这种差异，我想计算每一张图像的GLCM。

到目前为止，我使用的是"glcm“包，这是Alex Zvoleff开发的一个用于纹理分析的奇妙的包。不幸的是，它只能与滑动/移动窗口…一起使用但是因为我想为每个统计度量的每个图像接收一个值，所以对我来说似乎没有用……有谁可以帮助像我这样的R-新秀解决这个问题吗？:)

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

这个建议可能提供了通过EBImage包获得答案所需的工具。完整的答案可能需要对这里演示的纹理分析的结果应用额外的数据简化技术和统计分析。

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

对于EBImage，需要使用二进制掩码来定义对象，以便进行后续分析。在这种情况下，整个图像(数组)似乎用作分析对象，因此创建了覆盖整个图像的二进制掩模，然后对其进行修改以复制该示例。

# 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特征从灰度共生矩阵计算旋转不变纹理属性。参数haralick.scales用于指定纹理图案的预期重复比例。默认情况下，使用c(1, 2)每隔1和2个像素查找一次重复。这里我只是把它限制在1个像素。

我必须承认，我在没有完全理解它的情况下使用它。一个有用的资源可能是post by Earl Glynn。此外，Bioconductor上关于计算哈拉里克特征的question answered提供了很难找到的大量信息。

# 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

在这里，我使用热图来可视化和组织13个Haralick参数。该图非常清楚地显示了图像1和3与图像2非常相似，但也可以看到图像1和3之间的差异。

用于这个热图的矩阵，特别是如果它是从更多的图像生成的，可以通过主成分分析进行缩放和进一步分析，以识别相关图像。

  heatmap(m)

​

​

要了解有关EBImage的更多信息，请参阅the package vignette。

Answer 2

My R package GLCMTextures主要用于处理像glcm这样的空间栅格数据，但它也应该能够做到这一点。您必须分别为四个移位c(1，0)，c(1，1)，c(0，1)，c(-1，1)中的每一个分别列出GLCM，然后对每种类型的纹理度量进行平均，以获得方向不变的度量。

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