Speed of glcm_textures

[micha-silver]

I did a speed comparison between GLCM::glcm_textures() and the older package glcm::glcm(). It seems that glcm_textures is about 6 times slower (!). Is this expected? (I did the quantize_raster() step before starting microbenchmark)

Microbenchmark results:

expr       min        lq      mean    median        uq      max neval cld
 textures_old <- f1_glcm(red_raster)  25.28074  25.81914  26.23389  26.1424  26.70709  27.1866    10  a 
 textures_new <- f2_GLCMTextures(red_quantized) 142.20192 142.26140 142.73769   142.3798 142.76346 144.9905    10   b

Here's the code that I ran:

library(terra)
library(GLCMTextures)
library(glcm)
library(microbenchmark)

Data_dir <- "Namibia/Output/"
stk_file <- file.path(Data_dir, "stk_full.tif")
stk <- rast(stk_file)
# Clip to small area for testing
stk_crop <- crop(stk, ext(322400, 323000, 7885600, 7886200))
red <- stk_crop$red
red_raster <- raster::raster(red)

red_quantized <- quantize_raster(r = red, n_levels = 16,
                                 method = "equal prob")

grey_levels <- 16
wind <- c(5, 5)
shift <- list(c(0,1), c(1,1), c(1,0), c(-1,1))
f1_glcm = function(x) {textures_old <- glcm::glcm(x = x,
                               statistics=c('variance',
                                          'homogeneity',
                                          'contrast',
                                          'dissimilarity',
                                          'entropy'),
                               shift=shift,
                               n_grey = grey_levels, window=wind)

                      names(textures_old) <- c("variance", "homogeneity",
                                                "contrast", "dissimilarity", "entropy")
                      return(textures_old)
                      }
f2_GLCMTextures = function(x) {textures_new <- glcm_textures(r = x,
                                              w=wind,
                                              shift=shift,
                                              metrics = c("glcm_variance",
                                                          "glcm_homogeneity",
                                                          "glcm_contrast",
                                                          "glcm_dissimilarity",
                                                          "glcm_entropy"),
                                              quantization = "none",
                                              n_levels=grey_levels, na.rm = TRUE)

                      names(textures_new) <- c("variance", "homogeneity",
                                             "contrast","dissimilarity", "entropy")
                      return(textures_new)
                    }

microbenchmark(
               textures_old <- f1_glcm(red_raster),
               textures_new <- f2_GLCMTextures(red_quantized),
               times=10
              )

Thanks, Micha

[ailich]

Hi @micha-silver. Thanks for the benchmarks. I had never compared the speed of the two so I was unaware, though it is not that surprising as my C++ knowledge is limited, so the code likely has several areas where things could be optimized. Although GLCMTextures is newer, it was not created as a successor to glcm, but is just an alternative implementation. At this point speed is a secondary goal and while the C++ is likely sub-optimal, it drastically increases the speed compared to an R only implementation. The real strength of GLCMTextures in my mind is that it provides a way to calculate GLCM texture measures in a way that is transparent to the user as to exactly what is happening.

A lot of the parameters and style are similar to the glcm package because when I started to use GLCM calculations in my own research I found the glcm package syntax very clear compared to other alternatives, and I liked that it could explicitly handle NA values. That being said, I ran into issues where the correlation texture measure values were clearly being calculated incorrectly as they were beyond the range of -1 to 1 (azvoleff/glcm#6 ). Additionally, the glcm package compares their results against texture measures determined using ENVI which uses a different method for tabulating GLCM, which is different than the way it is described in the original literature. This modification makes it so that the GLCM is not symmetric across the diagonal. This has implications for if you need 4 or 8 shifts to get directionally invariant shifts, and it also means that you can get two different values for mean (µi and µj) and variance (σi and σj), whereas if the GLCM is symmetric the i and j versions are equivalent. I am unsure whether glcm uses the original method or the ENVI version of tabulating the GLCM as it compares to ENVI surfaces (8 directions and utilizes values outside the focal window) but states that you only need to use 4 directions to get directionally invariant textures which is consistent with the original literature. It has been a while, but I also tested the values of texture measures calculated by the glcm package other than correlation against what I would expect with the standard approach and got different results, and I also at one point got a trial license for ENVI and compared glcm and ENVI textures and got different results between the two software when using surfaces other than the test raster included in the glcm package. This confusion is why I developed the GLCMTextures package and the results should match those in the examples contained within Dr. Mryka Hall-Beyer's GLCM texture tutorial, which can be confirmed with make_glcm and glcm_metrics.

Here's some example code that shows the two do not produce equivalent results:

library(raster)
library(terra)
library(GLCMTextures)
library(glcm)

r<- raster(rast(volcano, extent= ext(2667400, 2667400 + ncol(volcano)*10, 6478700, 6478700 + nrow(volcano)*10), crs = "EPSG:27200"))

t1<- glcm(r, statistics=c("mean", "variance", "homogeneity", "contrast", "dissimilarity", "entropy", "second_moment", "correlation"), n_grey= 32, window =c(3,3), shift= c(1, 1), na_opt = "any") 

head(sort(unique(values(t1$glcm_correlation)), decreasing = TRUE)) #have infinite correlation and correlation of 5.8
#Inf 5.839971 1.000000 1.000000 1.000000 1.000000
head(sort(unique(values(t1$glcm_correlation)), decreasing = FALSE)) #have negative infinite correlation
# -Inf -0.5735393 -0.5000000 -0.5000000 -0.4588315 -0.4335550

t2<- glcm_textures(r, w = c(3,3), n_levels = 32, shift = c(1,1), metrics = c("glcm_mean", "glcm_variance", "glcm_homogeneity", "glcm_contrast", "glcm_dissimilarity", "glcm_entropy", "glcm_ASM", "glcm_correlation"), quantization = "equal range", na.rm=FALSE)
head(sort(unique(values(t2$glcm_correlation)), decreasing = TRUE)) #correlation max is 1
#1.0000000 0.9407407 0.8987342 0.8730159 0.8571429 0.8571429
head(sort(unique(values(t2$glcm_correlation)), decreasing = FALSE)) #correlation min is -1
# -1.0000000 -0.7142857 -0.6666667 -0.6000000 -0.5652174 -0.5555556

plot(t1)

plot(t2)

plot(t1-t2) #difference

[micha-silver]

Thanks for taking the time to compose the very detailed answer.

[ailich]

Also, @micha-silver, the documentation for the glcm package states that "to calculate GLCM textures over 'all directions' (in the terminology of commonly used remote sensing software), use: shift=list(c(0,1), c(1,1), c(1,0), c(1,-1)). This will calculate the average GLCM texture using shifts of 0 degrees, 45 degrees, 90 degrees, and 135 degrees." However, that only works if the GLCM is constructed in a symmetric way (that is why ENVI averages over 8 directions; see GLCMTextures ReadME for more details). For example, based on this a 45 degree shift and a 225 degree shift should produce equivalent results but in the glcm package they do not, while in GLCMTextures they do.

library(raster)
library(terra)
library(glcm)
library(GLCMTextures)

r<- raster(rast(volcano, extent= ext(2667400, 2667400 + ncol(volcano)*10, 6478700, 6478700 + nrow(volcano)*10), crs = "EPSG:27200"))

t1a<- glcm(r, statistics=c("mean", "variance", "homogeneity", "contrast", "dissimilarity", "entropy", "second_moment", "correlation"), n_grey= 32, window =c(3,3), shift= c(1, 1), na_opt = "any") #45 degree shift

t1b<- glcm(r, statistics=c("mean", "variance", "homogeneity", "contrast", "dissimilarity", "entropy", "second_moment", "correlation"), n_grey= 32, window =c(3,3), shift= c(-1, -1), na_opt = "any") #225 degree shift

plot(t1a-t1b) #ab diff

t2a<- glcm_textures(r, w = c(3,3), n_levels = 32, shift = c(1,1), metrics = c("glcm_mean", "glcm_variance", "glcm_homogeneity", "glcm_contrast", "glcm_dissimilarity", "glcm_entropy", "glcm_ASM", "glcm_correlation"), quantization = "equal range", na.rm=FALSE)
t2b<- glcm_textures(r, w = c(3,3), n_levels = 32, shift = c(-1,-1), metrics = c("glcm_mean", "glcm_variance", "glcm_homogeneity", "glcm_contrast", "glcm_dissimilarity", "glcm_entropy", "glcm_ASM", "glcm_correlation"), quantization = "equal range", na.rm=FALSE)

t2a-t2b
# names      : glcm_mean, glcm_variance, glcm_homogeneity, glcm_contrast, glcm_dissimilarity, glcm_entropy, glcm_ASM, glcm_correlation 
# min values :         0,             0,                0,             0,                  0,            0,        0,                0 
# max values :         0,             0,                0,             0,                  0,            0,        0,                0 

[rhijmans]

I think that there is some efficiency to be gained here ; especially if you want one or only a few metric(s) and you want to average them. Now you compute all metrics and then throw away the ones that are not requested.

It seems to me that all of that could be done in one call to the C++ code; something like this:

 out <- terra::focalCpp(r, w=w, fun = C_glcm_textures_helper, w2=w, n_levels= n_levels, shift = shifts, na_rm=na.rm, avg=avg_shifts, metrics=metrics)

(but it may not be worth the headache).

I think you can also replace this:

  if(avg_shifts){
    output<- terra::rast()
    for (j in 1:n_layers) {
      out_layer<- terra::mean(do.call(c, lapply(out_list, terra::subset,j)), na.rm=na.rm) #Average across all shifts
      output<- c(output, out_layer, warn=FALSE) #Create new stack of averaged values
    }} else{
      output<- out_list[[1]]
    }

with this:

 if(avg_shifts){
   output <- app(sds(out_list), "mean")
} else {
   output <- out_list[[1]]
}

(but whether that matters much, I do not know).

[rhijmans]

Another suggestion. Here, instead of:

    NumericVector curr_textures = C_glcm_metrics(curr_GLCM);
    out(i, 0) = curr_textures["glcm_contrast"];
    out(i, 1) = curr_textures["glcm_dissimilarity"];
    out(i, 2) = curr_textures["glcm_homogeneity"];
    out(i, 3) = curr_textures["glcm_ASM"];
    out(i, 4) = curr_textures["glcm_entropy"];
    out(i, 5) = curr_textures["glcm_mean"];
    out(i, 6) = curr_textures["glcm_variance"];
    out(i, 7) = curr_textures["glcm_correlation"];

Can't you do

    out(i, _) = C_glcm_metrics(curr_GLCM);

And generally avoid indexing by names.

[ailich]

Thanks for the suggestion @rhijmans. I made the change for out(i, _) = C_glcm_metrics(curr_GLCM) and using app(sds(out_list), "mean") as those are much cleaner than my previous code. Avoiding calculating all the metrics no matter what is also probably a good idea. I haven't implemented that yet, but will try to get to that in the not too distant future. Does indexing by name slow things down? I was trying to use name indices to prevent sloppy mistakes and make the code easier to follow, but perhaps that was misguided.

[rhijmans]

Preventing mistakes should always be the #1 priority (but speed is also important). Those changes avoid making copies of the data and that can help a lot. I do not know what the cost is of looking up an index by name (for each cell!). It would seem expensive, but you never know, as a lot depends on how that is implemented in Rcpp, and if the compiler can figure out that it only needs do the matching once. I suspect the compiler cannot figure that out in this case.

[ailich]

Thanks, I also realized I put, out(i, _) = curr_textures rather than out(i, _) = C_glcm_metrics(curr_GLCM) which may have still been copying the data, but I just pushed a new commit changing it to what you suggested.

[micha-silver]

@ailich Thanks to both of you for the efforts. Let me know when to do some additional tests on my end...

[ailich]

Some areas where modifications may help increase speed:

• i_mat and j_mat are defined within glcm_metrics and therefore reconstructed in each iteration even though they are the same each time.
• Pij is just an arma::mat copy of GLCM. Instead have C_make_glcm return an arma::mat so it doesn't need to be copied each time.
• curr_window (a matrix of focal values) copies values from xw (a vector containing the current focal values). Figure out how to directly fill in curr_window by row with values from x[Rcpp::Range(start,end)].
• Symmetric matrix is constructed by "double counting" . May be faster to just count one way and then add the transpose (not sure if this will be faster).

[ailich]

@micha-silver, I made some changes to improve speed (reduced duplication and only requested metrics are calculated). It is still slower than glcm but seems substantially faster than what it was before. Here's results from a slightly modified version of your code. Old is glcm, med is what GLCMTextures did at previous comit, and new is the current version of GLCMTextures. Also, I changed na_opt to "ignore" in glcm so that both had to remove NA values.

Unit: seconds
                                              expr       min       lq      mean    median        uq       max neval cld
 textures_old <- f1_glcm(red_raster)  29.39981  33.2366  34.03191  34.05833  35.35176  37.02302    10 a  
 textures_med <- f2_GLCMTextures(red_quantized) 565.40512 661.9819 668.54011 679.62102 694.94316 700.49068    10   c
 textures_new <- f3_GLCMTextures(red_quantized) 198.09181 209.4646 227.69237 234.78047 240.09895 249.20328    10  b 

[micha-silver]

Hi @ailich: I tried to duplicate your test on my end, and as far as speed, I am finding similar results as you. Thanks for your work on this.

glcm_textures is now much faster than before (but still about 4-5 times slower than the older package glcm. I assume we can attribute this to glcm not strictly following the literature.)

Here, using a subset of my study area, In the microbenchmark result below, AddTextures_old called glcm::glcm, and AddTextures called GLCMTextures::glcm_textures() with parallel parameters.

dim(WV_sharpened)
[1] 6000 4000    4
                                                expr       min        lq
AddTextures_old(WV_sharpened, save_tifs, grey_levels)  429.7173  429.7173 ...
AddTextures(WV_sharpened, save_tifs, grey_levels) 2883.3580 2883.3580 ...

Here are the results when rerunning using our whole study area

> dim(WV_sharpened)
[1] 16364 12509     4

                                                  expr      min       lq
AddTextures_old(WV_sharpened, save_tifs, grey_levels) 1949.862 1949.862 ...
AddTextures(WV_sharpened, save_tifs, grey_levels) 8320.902 8320.902 ...

However, when I use the whole study area, I am still seeing that strange splitup of the resulting textures: (On the smaller test area this does not happen.) The orange area below is our original satellite image, and the green is the GLCMTextures output.

` sessionInfo() R version 4.1.2 (2021-11-01) Platform: x86_64-pc-linux-gnu (64-bit) Running under: Debian GNU/Linux 11 (bullseye)

Matrix products: default BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.13.so

attached base packages: [1] parallel stats graphics grDevices utils datasets methods
[8] base

other attached packages: [1] microbenchmark_1.4.9 terra_1.5-19 fusionImage_0.0.1
[4] raster_3.5-15 sp_1.4-6 GLCMTextures_0.3.3
[7] supercells_0.6.5 corrplot_0.92 doParallel_1.0.16
[10] iterators_1.0.13 foreach_1.5.1 ranger_0.13.1
[13] caret_6.0-90 lattice_0.20-45 remotes_2.4.2
[16] exactextractr_0.7.2 glcm_1.6.5 RStoolbox_0.2.6
[19] sf_1.0-6 forcats_0.5.1 stringr_1.4.0
[22] dplyr_1.0.7 purrr_0.3.4 readr_2.1.2
[25] tidyr_1.1.4 tibble_3.1.6 ggplot2_3.3.5
[28] tidyverse_1.3.1

loaded via a namespace (and not attached): [1] TH.data_1.1-0 colorspace_2.0-2 ellipsis_0.3.2
[4] class_7.3-19 rgdal_1.5-28 fs_1.5.2
[7] rstudioapi_0.13 proxy_0.4-26 listenv_0.8.0
[10] mvtnorm_1.1-3 prodlim_2019.11.13 fansi_1.0.2
[13] lubridate_1.8.0 xml2_1.3.3 codetools_0.2-18
[16] splines_4.1.2 R.methodsS3_1.8.1 jsonlite_1.7.3
[19] pROC_1.18.0 broom_0.7.10 dbplyr_2.1.1
[22] R.oo_1.24.0 rgeos_0.5-9 compiler_4.1.2
[25] httr_1.4.2 backports_1.4.1 assertthat_0.2.1
[28] Matrix_1.3-4 cli_3.1.1 tools_4.1.2
[31] gtable_0.3.0 glue_1.6.1 reshape2_1.4.4
[34] Rcpp_1.0.8 cellranger_1.1.0 vctrs_0.3.8
[37] nlme_3.1-153 timeDate_3043.102 gower_0.2.2
[40] globals_0.14.0 rvest_1.0.2 lifecycle_1.0.1
[43] XML_3.99-0.8 future_1.23.0 zoo_1.8-9
[46] MASS_7.3-54 scales_1.1.1 ipred_0.9-12
[49] hms_1.1.1 sandwich_3.0-1 curl_4.3.2
[52] geosphere_1.5-14 rpart_4.1-15 stringi_1.7.6
[55] e1071_1.7-9 lava_1.6.10 rlang_1.0.1
[58] pkgconfig_2.0.3 recipes_0.1.17 tidyselect_1.1.1
[61] parallelly_1.30.0 plyr_1.8.6 magrittr_2.0.2
[64] R6_2.5.1 generics_0.1.2 multcomp_1.4-17
[67] DBI_1.1.2 pillar_1.7.0 haven_2.4.3
[70] withr_2.4.3 units_0.8-0 survival_3.2-13
[73] nnet_7.3-16 future.apply_1.8.1 modelr_0.1.8
[76] crayon_1.4.2 KernSmooth_2.23-20 utf8_1.2.2
[79] tzdb_0.2.0 grid_4.1.2 readxl_1.3.1
[82] data.table_1.14.2 ModelMetrics_1.2.2.2 reprex_2.0.1
[85] digest_0.6.29 classInt_0.4-3 R.utils_2.11.0
[88] stats4_4.1.2 munsell_0.5.0
`

Regards, Micha

[ailich]

@micha-silver thanks for confirming that it's also faster for you now. I may also be able to make a bit more improvements by having the function that tallies the counts return an arma::mat rather than an Rcpp::IntegerMatrix as that'd eliminate the need to convert the matrix from integer to a numeric by multiplying by 1.0 and rather than calculating the sum of the matrix through a loop , I could use the accu function to get the sum of all elements. In regards to the clipping issue, I'm not entirely sure why it is only occurring on Linux in this example; however, I was able to reproduce something similar on Windows and it seems related (rspatial/terra#519).

[ailich]

@micha-silver, would you be able to run comparisons between the current main branch and the arma_test branch which can be installed with remotes::install_github("ailich/GLCMTextures", ref="arma_test"). That arma_test branch uses Armadillo library matrices the whole way through which reduces the need to convert from Rcpp objects to armadillo objects, and allows for the transpose to be added rather than doing double counting. I think using the transpose may increase performance particularly in cases where the number of grey levels and/or the window size is larger. I ran some preliminary tests and I'm seeing comparable to increased performance on my end. Would you be able to run a few comparisons, and if it is indeed an improvement I'll push it to the main branch. Maybe use your example above but use one shift to save time (as that averaging among shifts is done outside my C++ functions) and try 16 grey levels with a 5x5 window, 16 grey levels with a 15x15 window, 32 grey levels with a 5x5 window, and 32 grey levels with a 15x15 window?

[micha-silver]

@ailich: For sure Currently I'm doing some tests on the server (128GB RAM) to complete the picture. Then I'll try the above.

[micha-silver]

Here are the timings when running on the server. The full raster result did not have any of the "splitting" of the original raster that I was seeing on my desktop machine (16GB vs 128GB on the server), so this might be a memory issue.

Here I again called both glcm::glcm and GLCMTextures::glcm_textures, but configured for only a single shift direction. Still the older glcm is more than 10x faster.

Full raster, single shift direction, two GLCM variables,

expr  min
AddTextures(red_terra, texture_file, save_tifs, grey_levels)   9657.7057
AddTextures_old(red_raster, texture_file_old, save_tifs, grey_levels)   729.9866

Clipped raster, single shift direction, two variables:

                                                                  expr      min
AddTextures(red_terra, texture_file, save_tifs, grey_levels) 2731.791
AddTextures_old(red_raster, texture_file_old, save_tifs, grey_levels)  159.703

[ailich]

@micha-silver I'm finding very similar performance between the main and arma_test branch with the main branch coming in slightly faster in this example. I ran the following code using each version and got a mean of about 18.5 -19 s.

library(terra)
library(microbenchmark)
library(GLCMTextures)

tmpFiles(current = TRUE, orphan = TRUE, old = TRUE, remove = TRUE)

r<- rast(nrow=1000, ncol=100, crs="EPSG:32601")
set.seed(5)
values(r)<- sample(0:31,size = ncell(r), replace = TRUE)

microbenchmark(out<- glcm_textures(r, w = c(15,15), n_levels = 32, shift = c(1,0), quantization = "none", na.rm = TRUE),
               times=10)

# main branch
# Unit: seconds
# min            lq    mean   median    uq      max      neval
# 18.5695   18.63823 18.7197 18.66345 18.79583 19.07438    10

# arma_test branch
# Unit: seconds      
# min         lq     mean     median      uq    max      neval
# 18.28288 18.56405 18.91579 18.99305 19.34786 19.45589    10