Analysis of Brodatz standard texture images: pixel-based vs. texture metrics

• Agustin.Lobo@ictja.csic.es

• 20190225, 20210314

• (Data and required functions)[]

• NOTE: Texture metrics are computed with Batch_GLCMHaralick.ijm in Fiji (Improvement: use package GLCMTextures instead https://github.com/ailich/GLCMTextures )

miR <- "C:/Users/1/ownCloudR/Rcommon/Rutils"
#require(rgdal)
#require(raster)
require(terra)
#require(RStoolbox)
require(tidyterra)
require(ggplot2)
require(MASS)
imadir <- "../Haralick"

1. Data input

Test images from http://sipi.usc.edu/database/database.php?volume=textures
“1.1.01.tiff” “1.2.03.tiff” “1.3.01.tiff” “1.5.07.tiff”
Four (4) 2000 x 200 subescenes were created for each using imagettes_log.R

imas <- list.files(path=imadir, glob2rx("*.tif"))
r <- rast(file.path(imadir,"ima1_1.tif"))
#ggR(r,stretch = "lin")
s1 <- rast(list.files(path=imadir,glob2rx("ima1*.tif"),full.names=TRUE))
s2 <- rast(list.files(path=imadir,glob2rx("ima2*.tif"),full.names=TRUE))
s3 <- rast(list.files(path=imadir,glob2rx("ima3*.tif"),full.names=TRUE))
s4 <- rast(list.files(path=imadir,glob2rx("ima4*.tif"),full.names=TRUE))
s <- c(s1,s2,s3,s4)
crs(s) <- NULL
#
#ggR(s, layer=1:16, geom_raster = TRUE) + 
#  scale_fill_gradientn(colours = grey.colors(10), guide=FALSE) +
#  facet_wrap(~layerName, ncol=4) + theme_void()

ggplot() +
  geom_spatraster(data=s) +
  scale_fill_gradient(low="grey5", high="grey98", guide=FALSE) +
  coord_sf() +
  facet_wrap(~lyr, ncol = 4)

2. Per-pixel classification

We first try with a conventional pixel-based classification:

#mostly hacked from http://remote-sensing.org/unsupervised-classification-with-r/
## returns the values of the raster dataset and write them in a matrix. 
v <- values(s)
v2 <- as.vector(v)
k <- kmeans(v2,4)
sclass <- s
values(sclass) <- k$cluster
sclass <- as.factor(sclass)

#ggR(sclass, layer=1:16, geom_raster = TRUE) + 
#  scale_fill_gradientn(colours = terrain.colors(4)) +
#  facet_wrap(~layerName, ncol=4) + theme_void()

##ggR(sclass, layer=1:16, geom_raster = TRUE) + 
##  scale_fill_brewer(palette="Set1") +
##  facet_wrap(~layerName, ncol=4) + theme_void()

##ggR(sclass[[1]], geom_raster = TRUE) + 
##  scale_fill_gradientn(colours = terrain.colors(4)) +
##  facet_wrap(~layerName, ncol=4) + theme_void()

ggplot() +
  geom_spatraster(data=sclass) +
  scale_color_gradientn(colours = terrain.colors(4)) +
  coord_sf() +
  facet_wrap(~lyr, ncol = 4)

3. Analysis of Texture metrics

Processed with Batch_GLCMHaralick.ijm in Fiji * https://imagej.nih.gov/ij/plugins/texture.html * https://imagej.nih.gov/ij/macros/Batch_GLCM_Measure.txt

GLCMHaralick.ijmbuilds a GLCM matrix by slicing a 3x3 window. One single GLCM, and its resulting set of texture metrics, is produced for each imagette.
The metrics are arranged in a table:

texturesDF <- read.csv("../Haralick/Haralick.csv",stringsAsFactors = FALSE)
texturesDF$type <- as.character(rep(1:4,rep(4,4)))
head(texturesDF)

NA

2.1 PCA

texturesPCA <- prcomp(data.matrix(texturesDF[,-c(1,7)]))
#plot(cumsum(texturesPCA$sdev)*100/sum(texturesPCA$sdev))

ggplot(data=data.frame(texturesDF[,c(1,7)],texturesPCA$x)) +
  geom_point(aes(x=PC1,y=PC2,color=type))

2.2 LDA

texturesLDA <- lda(texturesPCA$x[,1:3],grouping=texturesDF$type)
#plot(texturesLDA)
cl <- lda(texturesPCA$x[,1:3],grouping=texturesDF$type,CV=TRUE)$class
texturesLDApred <- predict(texturesLDA)
#knitr::kable(table(texturesDF$type,cl))
ggplot(data=data.frame(texturesDF[,c(1,7)],texturesLDApred$x,class=texturesLDApred$class)) +
  geom_text(aes(x=LD1,y=LD2,label=type,color=class),size=4)

table(texturesDF$type,cl)

   cl
    1 2 3 4
  1 4 0 0 0
  2 0 4 0 0
  3 0 0 4 0
  4 0 0 0 4

2.3 Hierarchical clustering

textureshc <- hclust(dist(texturesPCA$x[,1:3]),method="complete")
plot(textureshc,label=paste(texturesDF$Image, texturesLDApred$class),hang=-1)

texturesLDhc2 <- hclust(dist(texturesLDApred$x),method="complete")
plot(texturesLDhc2,label=paste(texturesDF$Image, texturesLDApred$class),hang=-1)

2.4 K-means

Original values:

#Original values:
texturesKM <- kmeans(texturesDF[,2:6],4)
table(texturesDF$type,texturesKM$cluster)

   
    1 2 3 4
  1 4 0 0 0
  2 1 0 3 0
  3 0 4 0 0
  4 0 2 0 2

#knitr::kable(table(texturesDF$type,texturesKM$cluster))

kmeans with of all PCs

#kmeans of all PCs
texturesKM <- kmeans(texturesPCA$x,4)
table(texturesDF$type,texturesKM$cluster)

   
    1 2 3 4
  1 0 0 4 0
  2 0 3 1 0
  3 4 0 0 0
  4 2 0 0 2

#knitr::kable(table(texturesDF$type,texturesKM$cluster))

kmeans with PC1 to PC3

#kmeans with PC1 to PC3
texturesKM <- kmeans(texturesPCA$x[,1:3],4)
table(texturesDF$type,texturesKM$cluster)

   
    1 2 3 4
  1 0 0 4 0
  2 0 3 1 0
  3 4 0 0 0
  4 2 0 0 2

#knitr::kable(table(texturesDF$type,texturesKM$cluster))

kmeans with all LDC

#kmeans with all LDC
texturesKM <- kmeans(texturesLDApred$x[,1:3],4)
table(texturesDF$type,texturesKM$cluster)

   
    1 2 3 4
  1 0 4 0 0
  2 4 0 0 0
  3 0 0 4 0
  4 0 0 0 4

#knitr::kable(table(texturesDF$type,texturesKM$cluster))

Original values, with initalization

#Original values:
texturesKM <- kmeans(texturesDF[,2:6],texturesDF[seq(1,16,by=4),2:6])
table(texturesDF$type,texturesKM$cluster)

   
    1 2 3 4
  1 4 0 0 0
  2 1 3 0 0
  3 0 0 4 0
  4 0 0 2 2

#knitr::kable(table(texturesDF$type,texturesKM$cluster))

kmeans with P1 to PC3, and with initialization

#kmeans with P1 to PC3, with initialization
texturesKM <- kmeans(texturesPCA$x[,1:3],texturesPCA$x[seq(1,16,by=4),1:3])
table(texturesDF$type,texturesKM$cluster)

   
    1 2 3 4
  1 4 0 0 0
  2 1 3 0 0
  3 0 0 4 0
  4 0 0 2 2

#knitr::kable(table(texturesDF$type,texturesKM$cluster))

kmeans with LDC1 to LDC3, and with initialization

#kmeans with LDC1 to LDC3, with initialization
texturesKM <- kmeans(texturesLDApred$x[,1:3],centers=texturesLDApred$x[seq(1,16,by=4),1:3])
table(texturesDF$type,texturesKM$cluster)

   
    1 2 3 4
  1 4 0 0 0
  2 0 4 0 0
  3 0 0 4 0
  4 0 0 0 4

#knitr::kable(table(texturesDF$type,texturesKM$cluster))