diff --git a/DESCRIPTION b/DESCRIPTION index d00be1e..054e4be 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: GLCMTextures Title: GLCM Textures of Raster Layers -Version: 0.5.3 +Version: 0.6 Authors@R: person(given = "Alexander", family = "Ilich", diff --git a/R/RcppExports.R b/R/RcppExports.R index 2cfd8d3..fa78c2e 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -1,23 +1,15 @@ # Generated by using Rcpp::compileAttributes() -> do not edit by hand # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 -C_make_glcm_counts <- function(x, n_levels, shift, na_rm) { - .Call(`_GLCMTextures_C_make_glcm_counts`, x, n_levels, shift, na_rm) +C_make_glcm <- function(x, n_levels, shift, na_rm, nrow, ncol, normalize) { + .Call(`_GLCMTextures_C_make_glcm`, x, n_levels, shift, na_rm, nrow, ncol, normalize) } -C_make_glcm <- function(x, n_levels, shift, na_rm) { - .Call(`_GLCMTextures_C_make_glcm`, x, n_levels, shift, na_rm) +C_glcm_metrics <- function(Pij, i_mat, j_mat, i_minus_j, i_minus_j2, i_minus_j_abs, n_levels, metric_indices, impute_corr) { + .Call(`_GLCMTextures_C_glcm_metrics`, Pij, i_mat, j_mat, i_minus_j, i_minus_j2, i_minus_j_abs, n_levels, metric_indices, impute_corr) } -C_GLSV <- function(Pij, n_levels) { - .Call(`_GLCMTextures_C_GLSV`, Pij, n_levels) -} - -C_glcm_metrics <- function(Pij, i_mat, j_mat, n_levels, metrics, impute_corr) { - .Call(`_GLCMTextures_C_glcm_metrics`, Pij, i_mat, j_mat, n_levels, metrics, impute_corr) -} - -C_glcm_textures_helper <- function(x, w2, n_levels, shift, metrics, na_rm, impute_corr, ni, nw) { - .Call(`_GLCMTextures_C_glcm_textures_helper`, x, w2, n_levels, shift, metrics, na_rm, impute_corr, ni, nw) +C_glcm_textures_helper <- function(x, w2, n_levels, shift_list, metric_indices, na_rm, impute_corr, ni, nw) { + .Call(`_GLCMTextures_C_glcm_textures_helper`, x, w2, n_levels, shift_list, metric_indices, na_rm, impute_corr, ni, nw) } diff --git a/R/glcm_metrics.R b/R/glcm_metrics.R index db5efd6..96d4082 100644 --- a/R/glcm_metrics.R +++ b/R/glcm_metrics.R @@ -22,16 +22,6 @@ glcm_metrics<-function(GLCM, metrics= c("glcm_contrast", "glcm_dissimilarity", " if (any(!(metrics %in% all_metrics))){ stop("Error: Invlaid metric. Valid metrics include 'glcm_contrast', 'glcm_dissimilarity', 'glcm_homogeneity', 'glcm_ASM', 'glcm_entropy', 'glcm_mean', 'glcm_variance', 'glcm_correlation'") } - needed_metrics<- metrics - if(("glcm_variance" %in% needed_metrics) & (!("glcm_mean" %in% needed_metrics))){ - needed_metrics<- c(needed_metrics, "glcm_mean") - } - if(("glcm_correlation" %in% needed_metrics) & (!("glcm_mean" %in% needed_metrics))){ - needed_metrics<- c(needed_metrics, "glcm_mean") - } - if(("glcm_correlation" %in% needed_metrics) & (!("glcm_variance" %in% needed_metrics))){ - needed_metrics<- c(needed_metrics, "glcm_variance") - } #Some metrics needed if(!is.list(GLCM)){GLCM<- list(GLCM)} out<- vector(mode="list", length = length(GLCM)) @@ -39,10 +29,16 @@ glcm_metrics<-function(GLCM, metrics= c("glcm_contrast", "glcm_dissimilarity", " for (i in 1:length(GLCM)) { i_mat<- matrix(data= 0:(nrow(GLCM[[i]])-1), nrow= nrow(GLCM[[i]]), ncol=ncol(GLCM[[i]]), byrow=FALSE) j_mat<- matrix(data= 0:(ncol(GLCM[[i]])-1), nrow= nrow(GLCM[[i]]), ncol=ncol(GLCM[[i]]), byrow=TRUE) + i_minus_j<- i_mat-j_mat + n_levels=nrow(GLCM[[i]]) - # k_vals<- seq((2*n_levels)-1)-1 - out[[i]]<- C_glcm_metrics(GLCM[[i]], i_mat = i_mat, j_mat = j_mat, n_levels=n_levels, metrics = needed_metrics, impute_corr= impute_corr) - out[[i]]<- out[[i]][names(out[[i]]) %in% metrics] + + metric_indices<- match(metrics, all_metrics) -1 + + out[[i]]<- C_glcm_metrics(GLCM[[i]], i_mat = i_mat, j_mat = j_mat, i_minus_j = i_minus_j, + i_minus_j2 = i_minus_j^2, i_minus_j_abs=abs(i_minus_j), + n_levels=n_levels, metric_indices = metric_indices, impute_corr=impute_corr) + names(out[[i]])<- metrics } if(length(out)==1){ out<- out[[1]] diff --git a/R/glcm_textures.R b/R/glcm_textures.R index ca73dfa..d64d270 100644 --- a/R/glcm_textures.R +++ b/R/glcm_textures.R @@ -81,17 +81,6 @@ glcm_textures<- function(r, w = c(3,3), n_levels, shift=list(c(1,0), c(1,1), c(0 warning("n_levels is > 32. This may be very computationally expensive.") } - needed_metrics<- metrics - if(("glcm_variance" %in% needed_metrics) & (!("glcm_mean" %in% needed_metrics))){ - needed_metrics<- c(needed_metrics, "glcm_mean") - } - if(("glcm_correlation" %in% needed_metrics) & (!("glcm_mean" %in% needed_metrics))){ - needed_metrics<- c(needed_metrics, "glcm_mean") - } - if(("glcm_correlation" %in% needed_metrics) & (!("glcm_variance" %in% needed_metrics))){ - needed_metrics<- c(needed_metrics, "glcm_variance") - } #Some metrics needed to calculate others - if(!is.null(quantization)){ warning("Use of 'quantization' is deprecated. Instead use 'quant_method'") if(is.null(quant_method)){ @@ -100,8 +89,9 @@ glcm_textures<- function(r, w = c(3,3), n_levels, shift=list(c(1,0), c(1,1), c(0 if(quantization == "equal range"){quant_method<- "range"} } } + if(is.null(quant_method)){stop("quant_method is NULL. Specify as 'none', 'range' or 'prob'")} - out_list<- vector(mode = "list", length=length(shift)) + if(quant_method!="none"){ r<- quantize_raster(r = r, n_levels = n_levels, quant_method = quant_method, min_val = min_val, max_val = max_val, maxcell=maxcell, wopt=wopt) } else if(!terra::is.int(r)){ @@ -110,19 +100,11 @@ glcm_textures<- function(r, w = c(3,3), n_levels, shift=list(c(1,0), c(1,1), c(0 if((unlist(terra::global(r, fun = max, na.rm=TRUE)) > (n_levels-1)) | (unlist(terra::global(r, fun = min, na.rm=TRUE)) < 0)){ stop("Error: raster must have values between 0 and n_levels-1")} - out_list<- vector(mode = "list", length=length(shift)) - for (k in 1:length(shift)) { - out_list[[k]]<- terra::focalCpp(r, w=w, fun = C_glcm_textures_helper, w2=w, n_levels= n_levels, shift = shift[[k]], metrics = needed_metrics, na_rm=na.rm, impute_corr = impute_corr, fillvalue=NA, wopt=wopt) - out_list[[k]]<- terra::subset(out_list[[k]], metrics, wopt=wopt) #Subset from needed to requested metrics - } - n_layers<- length(metrics) - if(length(shift) > 1){ - output<- terra::app(terra::sds(out_list), fun=mean, na.rm=na.rm, wopt=wopt) - } else{ - output<- out_list[[1]] - } - names(output)<- metrics #preserve names in case they were lost + metric_indices<- match(metrics, all_metrics)-1 + + output<- terra::focalCpp(r, w=w, fun = C_glcm_textures_helper, w2=w, n_levels= n_levels, shift = shift, metric_indices = metric_indices, na_rm=na.rm, impute_corr = impute_corr, fillvalue=NA, wopt=wopt) + names(output)<- metrics #Add in names if(include_scale){names(output)<- paste0(names(output), "_", w[1],"x", w[2])} #Add scale to layer names diff --git a/R/make_glcm.R b/R/make_glcm.R index 9a44bbf..896e730 100644 --- a/R/make_glcm.R +++ b/R/make_glcm.R @@ -19,9 +19,15 @@ #' @export make_glcm<- function(x, n_levels, shift, na.rm = FALSE, normalize=TRUE){ - if(class(x)[1]!="matrix"){ - x<- as.matrix(x, wide=TRUE) - } # Convert to matrix + nr=nrow(x) + nc=ncol(x) + + if(class(x)[1]=="matrix"){ + x<- as.vector(t(x)) + } else{ + x<- as.vector(t(x)) + } + if(isTRUE(any(x > (n_levels-1))) | isTRUE(any(x < 0))){ stop("Error: x must have values between 0 and n_levels-1") } @@ -29,14 +35,8 @@ make_glcm<- function(x, n_levels, shift, na.rm = FALSE, normalize=TRUE){ GLCM<- vector(mode="list", length = length(shift)) for (i in 1:length(shift)) { - if(normalize){ - GLCM[[i]]<- C_make_glcm(x=x, n_levels=n_levels, shift=shift[[i]], na_rm=na.rm) - } else{ - GLCM[[i]]<- C_make_glcm_counts(x=x, n_levels=n_levels, shift=shift[[i]], na_rm=na.rm) - } + GLCM[[i]]<- C_make_glcm(x=x, n_levels=n_levels, shift=shift[[i]], na_rm=na.rm, nrow=nr, ncol=nc,normalize=normalize) } if(length(GLCM)==1){GLCM<- GLCM[[1]]} return(GLCM) } - - diff --git a/README.md b/README.md index 23123d8..0dc1a86 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,7 @@ README ================ Alexander Ilich -January 02, 2025 +April 29, 2025 diff --git a/inst/testdata/txt_ep32_NA.RDS b/inst/testdata/txt_ep32_NA.RDS new file mode 100644 index 0000000..4378cf0 Binary files /dev/null and b/inst/testdata/txt_ep32_NA.RDS differ diff --git a/man/figures/README-textures1-1.png b/man/figures/README-textures1-1.png index 66a3179..7a245f6 100644 Binary files a/man/figures/README-textures1-1.png and b/man/figures/README-textures1-1.png differ diff --git a/man/figures/README-textures3-1.png b/man/figures/README-textures3-1.png index 1d77d24..d36431d 100644 Binary files a/man/figures/README-textures3-1.png and b/man/figures/README-textures3-1.png differ diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 370bf2a..2abfcb2 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -11,87 +11,65 @@ Rcpp::Rostream& Rcpp::Rcout = Rcpp::Rcpp_cout_get(); Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); #endif -// C_make_glcm_counts -IntegerMatrix C_make_glcm_counts(IntegerMatrix x, int n_levels, IntegerVector shift, bool na_rm); -RcppExport SEXP _GLCMTextures_C_make_glcm_counts(SEXP xSEXP, SEXP n_levelsSEXP, SEXP shiftSEXP, SEXP na_rmSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< IntegerMatrix >::type x(xSEXP); - Rcpp::traits::input_parameter< int >::type n_levels(n_levelsSEXP); - Rcpp::traits::input_parameter< IntegerVector >::type shift(shiftSEXP); - Rcpp::traits::input_parameter< bool >::type na_rm(na_rmSEXP); - rcpp_result_gen = Rcpp::wrap(C_make_glcm_counts(x, n_levels, shift, na_rm)); - return rcpp_result_gen; -END_RCPP -} // C_make_glcm -arma::mat C_make_glcm(IntegerMatrix x, int n_levels, IntegerVector shift, bool na_rm); -RcppExport SEXP _GLCMTextures_C_make_glcm(SEXP xSEXP, SEXP n_levelsSEXP, SEXP shiftSEXP, SEXP na_rmSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< IntegerMatrix >::type x(xSEXP); - Rcpp::traits::input_parameter< int >::type n_levels(n_levelsSEXP); - Rcpp::traits::input_parameter< IntegerVector >::type shift(shiftSEXP); - Rcpp::traits::input_parameter< bool >::type na_rm(na_rmSEXP); - rcpp_result_gen = Rcpp::wrap(C_make_glcm(x, n_levels, shift, na_rm)); - return rcpp_result_gen; -END_RCPP -} -// C_GLSV -NumericVector C_GLSV(arma::mat Pij, int n_levels); -RcppExport SEXP _GLCMTextures_C_GLSV(SEXP PijSEXP, SEXP n_levelsSEXP) { +arma::mat C_make_glcm(const IntegerVector& x, const int n_levels, const IntegerVector& shift, const bool na_rm, const int nrow, const int ncol, const bool normalize); +RcppExport SEXP _GLCMTextures_C_make_glcm(SEXP xSEXP, SEXP n_levelsSEXP, SEXP shiftSEXP, SEXP na_rmSEXP, SEXP nrowSEXP, SEXP ncolSEXP, SEXP normalizeSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< arma::mat >::type Pij(PijSEXP); - Rcpp::traits::input_parameter< int >::type n_levels(n_levelsSEXP); - rcpp_result_gen = Rcpp::wrap(C_GLSV(Pij, n_levels)); + Rcpp::traits::input_parameter< const IntegerVector& >::type x(xSEXP); + Rcpp::traits::input_parameter< const int >::type n_levels(n_levelsSEXP); + Rcpp::traits::input_parameter< const IntegerVector& >::type shift(shiftSEXP); + Rcpp::traits::input_parameter< const bool >::type na_rm(na_rmSEXP); + Rcpp::traits::input_parameter< const int >::type nrow(nrowSEXP); + Rcpp::traits::input_parameter< const int >::type ncol(ncolSEXP); + Rcpp::traits::input_parameter< const bool >::type normalize(normalizeSEXP); + rcpp_result_gen = Rcpp::wrap(C_make_glcm(x, n_levels, shift, na_rm, nrow, ncol, normalize)); return rcpp_result_gen; END_RCPP } // C_glcm_metrics -NumericVector C_glcm_metrics(arma::mat Pij, arma::mat i_mat, arma::mat j_mat, int n_levels, CharacterVector metrics, bool impute_corr); -RcppExport SEXP _GLCMTextures_C_glcm_metrics(SEXP PijSEXP, SEXP i_matSEXP, SEXP j_matSEXP, SEXP n_levelsSEXP, SEXP metricsSEXP, SEXP impute_corrSEXP) { +NumericVector C_glcm_metrics(const arma::mat& Pij, const arma::mat& i_mat, const arma::mat& j_mat, const arma::mat& i_minus_j, const arma::mat& i_minus_j2, const arma::mat& i_minus_j_abs, const int n_levels, const IntegerVector& metric_indices, const bool impute_corr); +RcppExport SEXP _GLCMTextures_C_glcm_metrics(SEXP PijSEXP, SEXP i_matSEXP, SEXP j_matSEXP, SEXP i_minus_jSEXP, SEXP i_minus_j2SEXP, SEXP i_minus_j_absSEXP, SEXP n_levelsSEXP, SEXP metric_indicesSEXP, SEXP impute_corrSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< arma::mat >::type Pij(PijSEXP); - Rcpp::traits::input_parameter< arma::mat >::type i_mat(i_matSEXP); - Rcpp::traits::input_parameter< arma::mat >::type j_mat(j_matSEXP); - Rcpp::traits::input_parameter< int >::type n_levels(n_levelsSEXP); - Rcpp::traits::input_parameter< CharacterVector >::type metrics(metricsSEXP); - Rcpp::traits::input_parameter< bool >::type impute_corr(impute_corrSEXP); - rcpp_result_gen = Rcpp::wrap(C_glcm_metrics(Pij, i_mat, j_mat, n_levels, metrics, impute_corr)); + Rcpp::traits::input_parameter< const arma::mat& >::type Pij(PijSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type i_mat(i_matSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type j_mat(j_matSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type i_minus_j(i_minus_jSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type i_minus_j2(i_minus_j2SEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type i_minus_j_abs(i_minus_j_absSEXP); + Rcpp::traits::input_parameter< const int >::type n_levels(n_levelsSEXP); + Rcpp::traits::input_parameter< const IntegerVector& >::type metric_indices(metric_indicesSEXP); + Rcpp::traits::input_parameter< const bool >::type impute_corr(impute_corrSEXP); + rcpp_result_gen = Rcpp::wrap(C_glcm_metrics(Pij, i_mat, j_mat, i_minus_j, i_minus_j2, i_minus_j_abs, n_levels, metric_indices, impute_corr)); return rcpp_result_gen; END_RCPP } // C_glcm_textures_helper -NumericMatrix C_glcm_textures_helper(IntegerVector x, IntegerVector w2, int n_levels, IntegerVector shift, CharacterVector metrics, bool na_rm, bool impute_corr, size_t ni, size_t nw); -RcppExport SEXP _GLCMTextures_C_glcm_textures_helper(SEXP xSEXP, SEXP w2SEXP, SEXP n_levelsSEXP, SEXP shiftSEXP, SEXP metricsSEXP, SEXP na_rmSEXP, SEXP impute_corrSEXP, SEXP niSEXP, SEXP nwSEXP) { +NumericMatrix C_glcm_textures_helper(const IntegerVector& x, const IntegerVector& w2, const int& n_levels, const List& shift_list, const IntegerVector& metric_indices, const bool na_rm, const bool impute_corr, size_t ni, size_t nw); +RcppExport SEXP _GLCMTextures_C_glcm_textures_helper(SEXP xSEXP, SEXP w2SEXP, SEXP n_levelsSEXP, SEXP shift_listSEXP, SEXP metric_indicesSEXP, SEXP na_rmSEXP, SEXP impute_corrSEXP, SEXP niSEXP, SEXP nwSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< IntegerVector >::type x(xSEXP); - Rcpp::traits::input_parameter< IntegerVector >::type w2(w2SEXP); - Rcpp::traits::input_parameter< int >::type n_levels(n_levelsSEXP); - Rcpp::traits::input_parameter< IntegerVector >::type shift(shiftSEXP); - Rcpp::traits::input_parameter< CharacterVector >::type metrics(metricsSEXP); - Rcpp::traits::input_parameter< bool >::type na_rm(na_rmSEXP); - Rcpp::traits::input_parameter< bool >::type impute_corr(impute_corrSEXP); + Rcpp::traits::input_parameter< const IntegerVector& >::type x(xSEXP); + Rcpp::traits::input_parameter< const IntegerVector& >::type w2(w2SEXP); + Rcpp::traits::input_parameter< const int& >::type n_levels(n_levelsSEXP); + Rcpp::traits::input_parameter< const List& >::type shift_list(shift_listSEXP); + Rcpp::traits::input_parameter< const IntegerVector& >::type metric_indices(metric_indicesSEXP); + Rcpp::traits::input_parameter< const bool >::type na_rm(na_rmSEXP); + Rcpp::traits::input_parameter< const bool >::type impute_corr(impute_corrSEXP); Rcpp::traits::input_parameter< size_t >::type ni(niSEXP); Rcpp::traits::input_parameter< size_t >::type nw(nwSEXP); - rcpp_result_gen = Rcpp::wrap(C_glcm_textures_helper(x, w2, n_levels, shift, metrics, na_rm, impute_corr, ni, nw)); + rcpp_result_gen = Rcpp::wrap(C_glcm_textures_helper(x, w2, n_levels, shift_list, metric_indices, na_rm, impute_corr, ni, nw)); return rcpp_result_gen; END_RCPP } static const R_CallMethodDef CallEntries[] = { - {"_GLCMTextures_C_make_glcm_counts", (DL_FUNC) &_GLCMTextures_C_make_glcm_counts, 4}, - {"_GLCMTextures_C_make_glcm", (DL_FUNC) &_GLCMTextures_C_make_glcm, 4}, - {"_GLCMTextures_C_GLSV", (DL_FUNC) &_GLCMTextures_C_GLSV, 2}, - {"_GLCMTextures_C_glcm_metrics", (DL_FUNC) &_GLCMTextures_C_glcm_metrics, 6}, + {"_GLCMTextures_C_make_glcm", (DL_FUNC) &_GLCMTextures_C_make_glcm, 7}, + {"_GLCMTextures_C_glcm_metrics", (DL_FUNC) &_GLCMTextures_C_glcm_metrics, 9}, {"_GLCMTextures_C_glcm_textures_helper", (DL_FUNC) &_GLCMTextures_C_glcm_textures_helper, 9}, {NULL, NULL, 0} }; diff --git a/src/glcm_cpp_functions.cpp b/src/glcm_cpp_functions.cpp index f13e89b..7f47661 100644 --- a/src/glcm_cpp_functions.cpp +++ b/src/glcm_cpp_functions.cpp @@ -1,151 +1,219 @@ -#include +#include // [[Rcpp::depends(RcppArmadillo)]] + using namespace Rcpp; using namespace arma; -//Make GLCM (non-normalized) +// Make GLCM // [[Rcpp::export]] -IntegerMatrix C_make_glcm_counts(IntegerMatrix x, int n_levels, IntegerVector shift, bool na_rm){ - IntegerMatrix GLCM(n_levels, n_levels);//initialize GLCM - int nr= x.nrow(); - int nc= x.ncol(); - - if((!na_rm) && (is_true(any(is_na(x))))){ - GLCM.fill(NA_INTEGER); - return GLCM; - } //Return window of NA's if any vals in window are NA - - if(na_rm && (is_true(all(is_na(x))))){ - GLCM.fill(NA_INTEGER); - return GLCM; - } //Return window of NA's if all values in window is NA - IntegerVector focalval_neighborval(2); - focalval_neighborval.names() = CharacterVector::create("focal_val", "neighborval"); - for(int i=0; i < nr; ++i){ - for(int j=0; j < nc; ++j){ - focalval_neighborval["focal_val"] = x(i,j); //focal val - IntegerVector neighbor_idx = {i-shift[1], j+shift[0]}; - if((neighbor_idx[0] < nr) && (neighbor_idx[1] < nc) && (neighbor_idx[0] >= 0) && (neighbor_idx[1] >= 0)){ - focalval_neighborval["neighborval"] = x(neighbor_idx[0], neighbor_idx[1]); //neighbor val - if(is_false(any(is_na(focalval_neighborval)))){ - GLCM(focalval_neighborval["focal_val"],focalval_neighborval["neighborval"])++; - GLCM(focalval_neighborval["neighborval"],focalval_neighborval["focal_val"])++; - } - }}} - return GLCM; -} +arma::mat C_make_glcm(const IntegerVector& x, + const int n_levels, + const IntegerVector& shift, + const bool na_rm, + const int nrow, + const int ncol, + const bool normalize) { -//Make GLCM (normalized) -// [[Rcpp::export]] -arma::mat C_make_glcm(IntegerMatrix x, int n_levels, IntegerVector shift, bool na_rm){ - IntegerMatrix GLCM = C_make_glcm_counts(x, n_levels, shift, na_rm); //tabulate counts - arma::mat GLCM_Norm(n_levels, n_levels); - if(is_true(any(is_na(GLCM)))){ - GLCM_Norm.fill(NA_REAL); - return GLCM_Norm; - } //If GLCM is NA return NA - - double total= sum(GLCM); - for(int m=0; m < GLCM.size(); ++m){ - GLCM_Norm[m] = (GLCM[m]*1.0)/total; - } //Normalize - return GLCM_Norm; -} + arma::mat GLCM(n_levels, n_levels, arma::fill::zeros); // single matrix from the start -//Gray Level Sum Vector -// [[Rcpp::export]] -NumericVector C_GLSV(arma::mat Pij, int n_levels) { - NumericVector k_prob(2*n_levels-1); //Initialize at zero - for (int i=0; i < n_levels; ++i) { - for (int j=0; j <= i; ++j) { - int k = i+j; - if(i==j){ - k_prob[k]= k_prob[k] + Pij(i,j); //add probability to corresponding spot in k_prob - } else{ - k_prob[k]= k_prob[k] + 2*Pij(i,j); //Since only looping through lower triangle have to count off diagonal elements twice + for (int i = 0; i < nrow; ++i) { + for (int j = 0; j < ncol; ++j) { + int idx = i * ncol + j; + int focal = x[idx]; + int ni = i - shift[1]; + int nj = j + shift[0]; + if (ni >= 0 && ni < nrow && nj >= 0 && nj < ncol) { + int nidx = ni * ncol + nj; + int neighbor = x[nidx]; + if (focal != NA_INTEGER && neighbor != NA_INTEGER) { + GLCM(focal, neighbor) += 1.0; + } else if (!na_rm) { + GLCM.fill(arma::datum::nan); + return GLCM; + } } } } - return k_prob; + + // Make symmetric + GLCM += GLCM.t(); + + // Normalize in place if requested + if (normalize) { + double total = arma::accu(GLCM); + if (total > 0.0) { + GLCM /= total; + } else { + GLCM.fill(arma::datum::nan); + } + } + + return GLCM; } -//Calculate Texture Metrics // [[Rcpp::export]] -NumericVector C_glcm_metrics(arma::mat Pij, arma::mat i_mat, arma::mat j_mat, int n_levels, CharacterVector metrics, bool impute_corr){ - - NumericVector textures = rep(NA_REAL,metrics.length()); - textures.names() = metrics; - if(!is_finite(Pij)){ - return textures;} - if(in(CharacterVector::create("glcm_contrast"), metrics)){ - textures["glcm_contrast"] = accu(Pij % pow(i_mat-j_mat,2)); //Contrast= sum(P_ij*(i-j)^2) - } - if(in(CharacterVector::create("glcm_dissimilarity"), metrics)){ - textures["glcm_dissimilarity"]= accu(Pij % abs(i_mat-j_mat)); //Dissimilarity= sum(P_ij*|i-j|) - } - if(in(CharacterVector::create("glcm_homogeneity"), metrics)){ - textures["glcm_homogeneity"]= accu(Pij/(1+(pow(i_mat-j_mat,2)))); //Homogeneity= sum(P_ij / (1+(i-j)^2)) - } - if(in(CharacterVector::create("glcm_ASM"), metrics)){ - textures["glcm_ASM"]= accu(pow(Pij,2)); //ASM= sum(P_ij^2) - } - if(in(CharacterVector::create("glcm_mean"), metrics)){ - textures["glcm_mean"]= accu(Pij % i_mat); //mean= sum(i*(P_ij)) - } - if(in(CharacterVector::create("glcm_variance"), metrics)){ - textures["glcm_variance"]= accu(Pij % pow(i_mat-textures["glcm_mean"],2)); //varaince= sum(P_ij*(i-u)^2 +NumericVector C_glcm_metrics(const arma::mat& Pij, + const arma::mat& i_mat, + const arma::mat& j_mat, + const arma::mat& i_minus_j, + const arma::mat& i_minus_j2, + const arma::mat& i_minus_j_abs, + const int n_levels, + const IntegerVector& metric_indices, + const bool impute_corr) { + + int n_metrics = metric_indices.length(); + NumericVector textures(n_metrics, NA_REAL); // Initialize output + + // Flags for what we need to compute + bool need_mean = false, need_var = false; + + for (int k = 0; k < n_metrics; ++k) { + int metric = metric_indices[k]; + if (metric == 5) need_mean = true; + if (metric == 6 || metric == 7) { + need_mean = true; + need_var = true; } - if(in(CharacterVector::create("glcm_correlation"), metrics)){ - if((textures["glcm_variance"]==0) && impute_corr){ - textures["glcm_correlation"]=0; //only way for all values to be the same in symmetric GLCM (limit approaches zero corr) - } else{ - textures["glcm_correlation"]= accu(Pij % (((i_mat-textures["glcm_mean"]) % (j_mat-textures["glcm_mean"]))/(textures["glcm_variance"]))); //Correlation= sum(P_ij*[((i-u)*(j-u))/(var)]) + } + + // Mean and variance (only if needed) + double glcm_mean = NA_REAL; + if (need_mean) glcm_mean = arma::accu(Pij % i_mat); + + double glcm_variance = NA_REAL; + if (need_var) glcm_variance = arma::accu(Pij % arma::square(i_mat - glcm_mean)); + + // Compute requested metrics + for (int m = 0; m < n_metrics; ++m) { + int metric = metric_indices[m]; + switch (metric) { + case 0: // Contrast + textures[m] = arma::accu(Pij % i_minus_j2); + break; + case 1: // Dissimilarity + textures[m] = arma::accu(Pij % i_minus_j_abs); + break; + case 2: // Homogeneity + textures[m] = arma::accu(Pij / (1.0 + i_minus_j2)); + break; + case 3: // ASM + textures[m] = arma::accu(arma::square(Pij)); + break; + case 4: { // Entropy (optimized) + double entropy = 0.0; + for (arma::uword i = 0; i < Pij.n_elem; ++i) { + double p = Pij[i]; + if (p > 0.0) { + entropy -= p * std::log(p); // avoids log(0) + } + } + textures[m] = entropy; + break; + } + case 5: // Mean + textures[m] = glcm_mean; + break; + case 6: // Variance + textures[m] = glcm_variance; + break; + case 7: { // Correlation (optimized) + if (glcm_variance == 0.0) { + textures[m] = impute_corr ? 0.0 : NA_REAL; + } else { + double corr_sum = 0.0; + for (arma::uword i = 0; i < Pij.n_elem; ++i) { + corr_sum += Pij[i] * (i_mat[i] - glcm_mean) * (j_mat[i] - glcm_mean); + } + textures[m] = corr_sum / glcm_variance; + } + break; } } - if(in(CharacterVector::create("glcm_entropy"), metrics)){ - arma::mat glcm_entropy_mat = Pij % ((-1) * log(Pij)); - glcm_entropy_mat.replace(datum::nan,0.0); - textures["glcm_entropy"]= accu(glcm_entropy_mat); //Entropy= sum(P_ij * (-ln(P_ij))) ; 0*ln(0)=0 - } - //if(in(CharacterVector::create("glcm_SA"), metrics)){ - //NumericVector k_prob = C_GLSV(Pij, n_levels); //Gray Level Sum Vector - //textures["glcm_SA"]= sum(k_prob*k_vals); //GLCM_SumAverage= sum(k*k_prob) - //} - return(textures); + } + + return textures; } -//GLCM across matrix using sliding window (terra) +// Focal function // [[Rcpp::export]] -NumericMatrix C_glcm_textures_helper(IntegerVector x, IntegerVector w2, int n_levels, IntegerVector shift, CharacterVector metrics, bool na_rm, bool impute_corr, size_t ni, size_t nw){ +NumericMatrix C_glcm_textures_helper(const IntegerVector& x, + const IntegerVector& w2, + const int& n_levels, + const List& shift_list, + const IntegerVector& metric_indices, + const bool na_rm, + const bool impute_corr, + size_t ni, + size_t nw) { - NumericMatrix out(ni, metrics.length()); + // Initialize output matrix + NumericMatrix out(ni, metric_indices.length()); out.fill(NA_REAL); - colnames(out)= metrics; - - arma::mat i_mat(n_levels,n_levels); - arma::mat j_mat(n_levels,n_levels); - // NumericVector k_vals((2*n_levels)-1); - // for (int k = 1; k < k_vals.length(); ++k) { - // k_vals[k] = k; - // } // All possible values of k, where k=i+j - for(int i=0; i 0) { + accum_metrics[m] /= valid_shifts[m]; + } else { + accum_metrics[m] = NA_REAL; + } + } + } else { + accum_metrics = accum_metrics / n_shifts; + } + + out(i, _) = accum_metrics; } - return(out); + + return out; } diff --git a/tests/testthat/test-glcm_textures.R b/tests/testthat/test-glcm_textures.R index e8c3eb4..89d1255 100644 --- a/tests/testthat/test-glcm_textures.R +++ b/tests/testthat/test-glcm_textures.R @@ -54,3 +54,15 @@ test_that("glcm_textures er works", { txt_er<- values(txt_er, mat=TRUE) expect_equal(txt_er, txt_er_extpected) }) + +test_that("glcm_textures NA handling and ordering of metrics works", { + txt_ep32_NA_expected <- readRDS(system.file("testdata", "txt_ep32_NA.RDS", package = "GLCMTextures")) + r<- rast(volcano, extent= ext(2667400, 2667400 + ncol(volcano)*10, 6478700, 6478700 + nrow(volcano)*10), crs = "EPSG:27200") #Use preloaded volcano dataset as a raster + set.seed(5) + r[sample(x = 1:ncell(r), size = 100)]<- NA + txt_ep32_NA<- glcm_textures(r, w = c(3,5), n_levels = 32, quant_method = "prob", shift=list(c(1, 0), c(1, 1), c(0, 1), c(-1, 1)), + metrics = rev(c("glcm_contrast", "glcm_dissimilarity", "glcm_homogeneity", "glcm_ASM", "glcm_entropy", "glcm_mean", "glcm_variance", "glcm_correlation")), + na.rm=TRUE, impute_corr = TRUE) + txt_ep32_NA<- values(txt_ep32_NA, mat=TRUE) + expect_equal(txt_ep32_NA, txt_ep32_NA_expected) +})