miR <- "C:/Users/1/ownCloudR/Rcommon/Rutils"
#require(rgdal)
#require(raster)
require(terra)
require(RStoolbox)
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)
ggR(s, layer=1:16, geom_raster = TRUE) + 
  scale_fill_gradientn(colours = grey.colors(10), guide=FALSE) +
  facet_wrap(~layerName, ncol=4) + theme_void()

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
ggR(sclass, layer=1:16, geom_raster = TRUE) + 
  scale_fill_gradientn(colours = terrain.colors(4)) +
  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()

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 2 0 2 0
  2 0 2 2 0
  3 0 0 0 4
  4 0 0 0 4
#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 2 0 0 2
  2 0 2 0 2
  3 0 0 4 0
  4 0 0 4 0
#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 0 0 0 4
  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 0 0 0 4
  3 4 0 0 0
  4 0 0 4 0
#knitr::kable(table(texturesDF$type,texturesKM$cluster))

Original values, with initalization

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

kmeans with P1 to PC3, and with initialization

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

kmeans with LDC1 to LDC3, and with initialization

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
LS0tCnRpdGxlOiAiQW5hbHlzaXMgb2YgQnJvZGF0eiBzdGFuZGFyZCB0ZXh0dXJlIGltYWdlczogcGl4ZWwtYmFzZWQgdnMuIHRleHR1cmUgbWV0cmljcyAiCm91dHB1dDoKICBodG1sX25vdGVib29rOgogICAgY29kZV9mb2xkaW5nOiBoaWRlCiAgICBmaWdfY2FwdGlvbjogVFJVRQotLS0KCiogQWd1c3Rpbi5Mb2JvQGljdGphLmNzaWMuZXMKKiAyMDE5MDIyNSwgMjAyMTAzMTQKKiAoRGF0YSBhbmQgcmVxdWlyZWQgZnVuY3Rpb25zKVtdIAoKYGBge3IgfQptaVIgPC0gIkM6L1VzZXJzLzEvb3duQ2xvdWRSL1Jjb21tb24vUnV0aWxzIgojcmVxdWlyZShyZ2RhbCkKI3JlcXVpcmUocmFzdGVyKQpyZXF1aXJlKHRlcnJhKQpyZXF1aXJlKFJTdG9vbGJveCkKcmVxdWlyZShnZ3Bsb3QyKQpyZXF1aXJlKE1BU1MpCmltYWRpciA8LSAiLi4vSGFyYWxpY2siCmBgYAoKIyMgMS4gRGF0YSBpbnB1dAoKVGVzdCBpbWFnZXMgZnJvbSBodHRwOi8vc2lwaS51c2MuZWR1L2RhdGFiYXNlL2RhdGFiYXNlLnBocD92b2x1bWU9dGV4dHVyZXMgIAoiMS4xLjAxLnRpZmYiICIxLjIuMDMudGlmZiIgIjEuMy4wMS50aWZmIiAiMS41LjA3LnRpZmYiICAgCkZvdXIgKDQpIDIwMDAgeCAyMDAgc3ViZXNjZW5lcyB3ZXJlIGNyZWF0ZWQgZm9yIGVhY2ggdXNpbmcgaW1hZ2V0dGVzX2xvZy5SCgpgYGB7ciBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFLCByZXN1bHRzPSJob2xkIiwgZmlnLndpZHRoPTgsIGZpZy5oZWlnaHQ9MTB9CmltYXMgPC0gbGlzdC5maWxlcyhwYXRoPWltYWRpciwgZ2xvYjJyeCgiKi50aWYiKSkKciA8LSByYXN0KGZpbGUucGF0aChpbWFkaXIsImltYTFfMS50aWYiKSkKI2dnUihyLHN0cmV0Y2ggPSAibGluIikKczEgPC0gcmFzdChsaXN0LmZpbGVzKHBhdGg9aW1hZGlyLGdsb2IycngoImltYTEqLnRpZiIpLGZ1bGwubmFtZXM9VFJVRSkpCnMyIDwtIHJhc3QobGlzdC5maWxlcyhwYXRoPWltYWRpcixnbG9iMnJ4KCJpbWEyKi50aWYiKSxmdWxsLm5hbWVzPVRSVUUpKQpzMyA8LSByYXN0KGxpc3QuZmlsZXMocGF0aD1pbWFkaXIsZ2xvYjJyeCgiaW1hMyoudGlmIiksZnVsbC5uYW1lcz1UUlVFKSkKczQgPC0gcmFzdChsaXN0LmZpbGVzKHBhdGg9aW1hZGlyLGdsb2IycngoImltYTQqLnRpZiIpLGZ1bGwubmFtZXM9VFJVRSkpCnMgPC0gYyhzMSxzMixzMyxzNCkKZ2dSKHMsIGxheWVyPTE6MTYsIGdlb21fcmFzdGVyID0gVFJVRSkgKyAKICBzY2FsZV9maWxsX2dyYWRpZW50bihjb2xvdXJzID0gZ3JleS5jb2xvcnMoMTApLCBndWlkZT1GQUxTRSkgKwogIGZhY2V0X3dyYXAofmxheWVyTmFtZSwgbmNvbD00KSArIHRoZW1lX3ZvaWQoKQpgYGAKIyMgMi4gUGVyLXBpeGVsIGNsYXNzaWZpY2F0aW9uCgpXZSBmaXJzdCB0cnkgd2l0aCBhIGNvbnZlbnRpb25hbCBwaXhlbC1iYXNlZCBjbGFzc2lmaWNhdGlvbjoKCmBgYHtyIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0UsIHJlc3VsdHM9ImhvbGQiLCBmaWcud2lkdGg9OCwgZmlnLmhlaWdodD0xMH0KI21vc3RseSBoYWNrZWQgZnJvbSBodHRwOi8vcmVtb3RlLXNlbnNpbmcub3JnL3Vuc3VwZXJ2aXNlZC1jbGFzc2lmaWNhdGlvbi13aXRoLXIvCiMjIHJldHVybnMgdGhlIHZhbHVlcyBvZiB0aGUgcmFzdGVyIGRhdGFzZXQgYW5kIHdyaXRlIHRoZW0gaW4gYSBtYXRyaXguIAp2IDwtIHZhbHVlcyhzKQp2MiA8LSBhcy52ZWN0b3IodikKayA8LSBrbWVhbnModjIsNCkKc2NsYXNzIDwtIHMKdmFsdWVzKHNjbGFzcykgPC0gayRjbHVzdGVyCmdnUihzY2xhc3MsIGxheWVyPTE6MTYsIGdlb21fcmFzdGVyID0gVFJVRSkgKyAKICBzY2FsZV9maWxsX2dyYWRpZW50bihjb2xvdXJzID0gdGVycmFpbi5jb2xvcnMoNCkpICsKICBmYWNldF93cmFwKH5sYXllck5hbWUsIG5jb2w9NCkgKyB0aGVtZV92b2lkKCkKCiNnZ1Ioc2NsYXNzW1sxXV0sIGdlb21fcmFzdGVyID0gVFJVRSkgKyAKIyAgc2NhbGVfZmlsbF9ncmFkaWVudG4oY29sb3VycyA9IHRlcnJhaW4uY29sb3JzKDQpKSArCiMgIGZhY2V0X3dyYXAofmxheWVyTmFtZSwgbmNvbD00KSArIHRoZW1lX3ZvaWQoKQpgYGAKCgojIyAzLiBBbmFseXNpcyBvZiBUZXh0dXJlIG1ldHJpY3MKUHJvY2Vzc2VkIHdpdGggQmF0Y2hfR0xDTUhhcmFsaWNrLmlqbSBpbiBGaWppCiogaHR0cHM6Ly9pbWFnZWoubmloLmdvdi9pai9wbHVnaW5zL3RleHR1cmUuaHRtbAoqIGh0dHBzOi8vaW1hZ2VqLm5paC5nb3YvaWovbWFjcm9zL0JhdGNoX0dMQ01fTWVhc3VyZS50eHQKCkdMQ01IYXJhbGljay5pam1idWlsZHMgYSBHTENNIG1hdHJpeCBieSBzbGljaW5nIGEgM3gzIHdpbmRvdy4gT25lIHNpbmdsZSBHTENNLCBhbmQgaXRzIHJlc3VsdGluZyBzZXQgb2YgdGV4dHVyZSBtZXRyaWNzLCBpcyBwcm9kdWNlZCBmb3IgZWFjaCBpbWFnZXR0ZS4gIApUaGUgbWV0cmljcyBhcmUgYXJyYW5nZWQgaW4gYSB0YWJsZToKCmBgYHtyIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0V9CnRleHR1cmVzREYgPC0gcmVhZC5jc3YoIi4uL0hhcmFsaWNrL0hhcmFsaWNrLmNzdiIsc3RyaW5nc0FzRmFjdG9ycyA9IEZBTFNFKQp0ZXh0dXJlc0RGJHR5cGUgPC0gYXMuY2hhcmFjdGVyKHJlcCgxOjQscmVwKDQsNCkpKQpoZWFkKHRleHR1cmVzREYpCgpgYGAKCiMjIyAyLjEgUENBCgpgYGB7ciBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFfQp0ZXh0dXJlc1BDQSA8LSBwcmNvbXAoZGF0YS5tYXRyaXgodGV4dHVyZXNERlssLWMoMSw3KV0pKQojcGxvdChjdW1zdW0odGV4dHVyZXNQQ0Ekc2RldikqMTAwL3N1bSh0ZXh0dXJlc1BDQSRzZGV2KSkKCmdncGxvdChkYXRhPWRhdGEuZnJhbWUodGV4dHVyZXNERlssYygxLDcpXSx0ZXh0dXJlc1BDQSR4KSkgKwogIGdlb21fcG9pbnQoYWVzKHg9UEMxLHk9UEMyLGNvbG9yPXR5cGUpKQpgYGAKCiMjIyMgMi4yIExEQQoKYGBge3IgbWVzc2FnZT1GQUxTRSwgd2FybmluZz1GQUxTRX0KdGV4dHVyZXNMREEgPC0gbGRhKHRleHR1cmVzUENBJHhbLDE6M10sZ3JvdXBpbmc9dGV4dHVyZXNERiR0eXBlKQojcGxvdCh0ZXh0dXJlc0xEQSkKY2wgPC0gbGRhKHRleHR1cmVzUENBJHhbLDE6M10sZ3JvdXBpbmc9dGV4dHVyZXNERiR0eXBlLENWPVRSVUUpJGNsYXNzCnRleHR1cmVzTERBcHJlZCA8LSBwcmVkaWN0KHRleHR1cmVzTERBKQoja25pdHI6OmthYmxlKHRhYmxlKHRleHR1cmVzREYkdHlwZSxjbCkpCmdncGxvdChkYXRhPWRhdGEuZnJhbWUodGV4dHVyZXNERlssYygxLDcpXSx0ZXh0dXJlc0xEQXByZWQkeCxjbGFzcz10ZXh0dXJlc0xEQXByZWQkY2xhc3MpKSArCiAgZ2VvbV90ZXh0KGFlcyh4PUxEMSx5PUxEMixsYWJlbD10eXBlLGNvbG9yPWNsYXNzKSxzaXplPTQpCnRhYmxlKHRleHR1cmVzREYkdHlwZSxjbCkKYGBgCgojIyMgMi4zIEhpZXJhcmNoaWNhbCBjbHVzdGVyaW5nCgpgYGB7ciB9CnRleHR1cmVzaGMgPC0gaGNsdXN0KGRpc3QodGV4dHVyZXNQQ0EkeFssMTozXSksbWV0aG9kPSJjb21wbGV0ZSIpCnBsb3QodGV4dHVyZXNoYyxsYWJlbD1wYXN0ZSh0ZXh0dXJlc0RGJEltYWdlLCB0ZXh0dXJlc0xEQXByZWQkY2xhc3MpLGhhbmc9LTEpCnRleHR1cmVzTERoYzIgPC0gaGNsdXN0KGRpc3QodGV4dHVyZXNMREFwcmVkJHgpLG1ldGhvZD0iY29tcGxldGUiKQpwbG90KHRleHR1cmVzTERoYzIsbGFiZWw9cGFzdGUodGV4dHVyZXNERiRJbWFnZSwgdGV4dHVyZXNMREFwcmVkJGNsYXNzKSxoYW5nPS0xKQpgYGAKCiMjIyAyLjQgSy1tZWFucwpPcmlnaW5hbCB2YWx1ZXM6CmBgYHtyLCByZXN1bHRzID0gImhvbGQifQojT3JpZ2luYWwgdmFsdWVzOgp0ZXh0dXJlc0tNIDwtIGttZWFucyh0ZXh0dXJlc0RGWywyOjZdLDQpCnRhYmxlKHRleHR1cmVzREYkdHlwZSx0ZXh0dXJlc0tNJGNsdXN0ZXIpCiNrbml0cjo6a2FibGUodGFibGUodGV4dHVyZXNERiR0eXBlLHRleHR1cmVzS00kY2x1c3RlcikpCmBgYAprbWVhbnMgd2l0aCBvZiBhbGwgUENzCmBgYHtyICwgcmVzdWx0cyA9ICJob2xkIn0KI2ttZWFucyBvZiBhbGwgUENzCnRleHR1cmVzS00gPC0ga21lYW5zKHRleHR1cmVzUENBJHgsNCkKdGFibGUodGV4dHVyZXNERiR0eXBlLHRleHR1cmVzS00kY2x1c3RlcikKI2tuaXRyOjprYWJsZSh0YWJsZSh0ZXh0dXJlc0RGJHR5cGUsdGV4dHVyZXNLTSRjbHVzdGVyKSkKYGBgCmttZWFucyB3aXRoIFBDMSB0byBQQzMKYGBge3IgLCByZXN1bHRzID0gImhvbGQifQoja21lYW5zIHdpdGggUEMxIHRvIFBDMwp0ZXh0dXJlc0tNIDwtIGttZWFucyh0ZXh0dXJlc1BDQSR4WywxOjNdLDQpCnRhYmxlKHRleHR1cmVzREYkdHlwZSx0ZXh0dXJlc0tNJGNsdXN0ZXIpCiNrbml0cjo6a2FibGUodGFibGUodGV4dHVyZXNERiR0eXBlLHRleHR1cmVzS00kY2x1c3RlcikpCmBgYAprbWVhbnMgd2l0aCBhbGwgTERDCmBgYHtyICwgcmVzdWx0cyA9ICJob2xkIn0KI2ttZWFucyB3aXRoIGFsbCBMREMKdGV4dHVyZXNLTSA8LSBrbWVhbnModGV4dHVyZXNMREFwcmVkJHhbLDE6M10sNCkKdGFibGUodGV4dHVyZXNERiR0eXBlLHRleHR1cmVzS00kY2x1c3RlcikKI2tuaXRyOjprYWJsZSh0YWJsZSh0ZXh0dXJlc0RGJHR5cGUsdGV4dHVyZXNLTSRjbHVzdGVyKSkKYGBgCgpPcmlnaW5hbCB2YWx1ZXMsIHdpdGggaW5pdGFsaXphdGlvbgpgYGB7ciwgcmVzdWx0cyA9ICJob2xkIn0KI09yaWdpbmFsIHZhbHVlczoKdGV4dHVyZXNLTSA8LSBrbWVhbnModGV4dHVyZXNERlssMjo2XSx0ZXh0dXJlc0RGW3NlcSgxLDE2LGJ5PTQpLDI6Nl0pCnRhYmxlKHRleHR1cmVzREYkdHlwZSx0ZXh0dXJlc0tNJGNsdXN0ZXIpCiNrbml0cjo6a2FibGUodGFibGUodGV4dHVyZXNERiR0eXBlLHRleHR1cmVzS00kY2x1c3RlcikpCmBgYAoKa21lYW5zIHdpdGggUDEgdG8gUEMzLCBhbmQgd2l0aCBpbml0aWFsaXphdGlvbgpgYGB7ciwgcmVzdWx0cyA9ICJob2xkIiB9CiNrbWVhbnMgd2l0aCBQMSB0byBQQzMsIHdpdGggaW5pdGlhbGl6YXRpb24KdGV4dHVyZXNLTSA8LSBrbWVhbnModGV4dHVyZXNQQ0EkeFssMTozXSx0ZXh0dXJlc1BDQSR4W3NlcSgxLDE2LGJ5PTQpLDE6M10pCnRhYmxlKHRleHR1cmVzREYkdHlwZSx0ZXh0dXJlc0tNJGNsdXN0ZXIpCiNrbml0cjo6a2FibGUodGFibGUodGV4dHVyZXNERiR0eXBlLHRleHR1cmVzS00kY2x1c3RlcikpCmBgYAprbWVhbnMgd2l0aCBMREMxIHRvIExEQzMsIGFuZCB3aXRoIGluaXRpYWxpemF0aW9uCmBgYHtyICwgcmVzdWx0cyA9ICJob2xkIn0KI2ttZWFucyB3aXRoIExEQzEgdG8gTERDMywgd2l0aCBpbml0aWFsaXphdGlvbgp0ZXh0dXJlc0tNIDwtIGttZWFucyh0ZXh0dXJlc0xEQXByZWQkeFssMTozXSxjZW50ZXJzPXRleHR1cmVzTERBcHJlZCR4W3NlcSgxLDE2LGJ5PTQpLDE6M10pCnRhYmxlKHRleHR1cmVzREYkdHlwZSx0ZXh0dXJlc0tNJGNsdXN0ZXIpCiNrbml0cjo6a2FibGUodGFibGUodGV4dHVyZXNERiR0eXBlLHRleHR1cmVzS00kY2x1c3RlcikpCmBgYAoKCg==