diff --git a/DESCRIPTION b/DESCRIPTION index fab7c9a..023d39b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: RStoolbox Type: Package Title: Remote Sensing Data Analysis -Version: 1.0.0 -Date: 2024-04-24 +Version: 1.0.1 +Date: 2024-08-26 Authors@R: c( person("Benjamin", "Leutner", role= "aut", email="rstoolboxpackage@gmail.com", comment = c(ORCID = "0000-0002-6893-2002")), person("Ned", "Horning", role ="aut", email="horning@amnh.org"), @@ -34,6 +34,7 @@ Imports: magrittr Suggests: randomForest, + lattice, kernlab, e1071, gridExtra, @@ -42,5 +43,5 @@ Suggests: LinkingTo: Rcpp, RcppArmadillo License: GPL (>=3) -RoxygenNote: 7.3.1 +RoxygenNote: 7.3.2 LazyData: true diff --git a/R/RcppExports.R b/R/RcppExports.R index 80758d0..8c277fc 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -49,7 +49,7 @@ specSimC <- function(x, em) { .Call('_RStoolbox_specSimC', PACKAGE = 'RStoolbox', x, em) } -spectralIndicesCpp <- function(x, indices, redBand, blueBand, greenBand, nirBand, redEdge1Band, redEdge2Band, redEdge3Band, swir1Band, swir2Band, swir3Band, maskLayer, maskValue, L, s, G, C1, C2, Levi, swir2ccc, swir2cdiff, sf) { - .Call('_RStoolbox_spectralIndicesCpp', PACKAGE = 'RStoolbox', x, indices, redBand, blueBand, greenBand, nirBand, redEdge1Band, redEdge2Band, redEdge3Band, swir1Band, swir2Band, swir3Band, maskLayer, maskValue, L, s, G, C1, C2, Levi, swir2ccc, swir2cdiff, sf) +spectralIndicesCpp <- function(x, indices, redBand, blueBand, greenBand, nirBand, redEdge1Band, redEdge2Band, redEdge3Band, swir1Band, swir2Band, swir3Band, maskLayer, maskValue, L, s, G, C1, C2, Levi, swir2ccc, swir2cdiff, sf, formulas) { + .Call('_RStoolbox_spectralIndicesCpp', PACKAGE = 'RStoolbox', x, indices, redBand, blueBand, greenBand, nirBand, redEdge1Band, redEdge2Band, redEdge3Band, swir1Band, swir2Band, swir3Band, maskLayer, maskValue, L, s, G, C1, C2, Levi, swir2ccc, swir2cdiff, sf, formulas) } diff --git a/R/internalFunctions.R b/R/internalFunctions.R index 611eb77..82ea28d 100644 --- a/R/internalFunctions.R +++ b/R/internalFunctions.R @@ -255,7 +255,84 @@ #' On package startup #' @noRd .onLoad <- function(libname, pkgname){ + + L <- 0.5 + G <- 2.5 + L_evi <- 1 + C1 <- 6 + C2 <- 7.5 + s <- 1 + swir2ccc <- NULL + swir2coc <- NULL + NDVI <- NULL + + .IDXDB <- list( + CLG = list(c("Gitelson2003", "Green-band Chlorophyll Index"), + function(redEdge3, green) {redEdge3/green - 1}), + CLRE = list(c("Gitelson2003", "Red-edge-band Chlorophyll Index"), + function(redEdge3, redEdge1) {redEdge3/redEdge1 - 1}), + CTVI = list(c("Perry1984", "Corrected Transformed Vegetation Index"), + function(red, nir) {(NDVI+.5)/sqrt(abs(NDVI+.5))} ), + DVI = list(c("Richardson1977", "Difference Vegetation Index"), + function(red, nir) {s*nir-red}), + EVI = list(c("Huete1999", "Enhanced Vegetation Index"), + function(red, nir, blue) {G * ((nir - red) / (nir + C1 * red - C2 * blue + L_evi))}), + EVI2 = list(c("Jiang 2008", "Two-band Enhanced Vegetation Index"), + function(red, nir) {G * (nir-red)/(nir + 2.4*red +1)}), + GEMI = list(c("Pinty1992", "Global Environmental Monitoring Index"), + function(red, nir) {(((nir^2 - red^2) * 2 + (nir * 1.5) + (red * 0.5) ) / (nir + red + 0.5)) * (1 - ((((nir^2 - red^2) * 2 + (nir * 1.5) + (red * 0.5) ) / (nir + red + 0.5)) * 0.25)) - ((red - 0.125) / (1 - red))}), + GNDVI = list(c("Gitelson1998", "Green Normalised Difference Vegetation Index" ), + function(green, nir) {(nir-green)/(nir+green)}), + KNDVI = list(c("Camps-Valls2021", "Kernel Normalised Difference Vegetation Index"), + function(red, nir) {tanh(((nir-red)/(nir+red)))^2}), + MCARI = list(c("Daughtery2000", "Modified Chlorophyll Absorption Ratio Index" ), + function(green, red, redEdge1) {((redEdge1-red)-(redEdge1-green))*(redEdge1/red)}), + MNDWI = list(c("Xu2006", "Modified Normalised Difference Water Index"), + function(green, swir2) {(green-swir2) / (green+swir2)}), + MSAVI = list(c("Qi1994", "Modified Soil Adjusted Vegetation Index" ), + function(red, nir) {nir + 0.5 - (0.5 * sqrt((2 * nir + 1)^2 - 8 * (nir - (2 * red))))}), + MSAVI2 = list(c("Qi1994", "Modified Soil Adjusted Vegetation Index 2" ), + function(red, nir) {(2 * (nir + 1) - sqrt((2 * nir + 1)^2 - 8 * (nir - red))) / 2}), + MTCI = list(c("DashAndCurran2004", "MERIS Terrestrial Chlorophyll Index"), + function(red, redEdge1, redEdge2) {(redEdge2-redEdge1)/(redEdge1-red)}), + NBRI = list(c("Garcia1991", "Normalised Burn Ratio Index"), + function(nir, swir3) { (nir - swir3) / (nir + swir3)}), + NDREI1 = list(c("GitelsonAndMerzlyak1994", "Normalised Difference Red Edge Index 1"), + function(redEdge2, redEdge1) {(redEdge2-redEdge1)/(redEdge2+redEdge1)}), + NDREI2 = list(c("Barnes2000", "Normalised Difference Red Edge Index 2"), + function(redEdge3, redEdge1) {(redEdge3-redEdge1)/(redEdge3+redEdge1)}), + NDVI = list(c("Rouse1974", "Normalised Difference Vegetation Index"), + function(red, nir) {(nir-red)/(nir+red)}), + NDVIC = list(c("Nemani1993", "Corrected Normalised Difference Vegetation Index" ), + function(red, nir, swir2) {(nir-red)/(nir+red)*(1-((swir2 - swir2ccc)/(swir2coc-swir2ccc)))}), + NDWI = list(c("McFeeters1996", "Normalised Difference Water Index"), + function(green, nir) {(green - nir)/(green + nir)}), + NDWI2 = list(c("Gao1996", "Normalised Difference Water Index"), + function(nir, swir2) {(nir - swir2)/(nir + swir2)}), + NRVI = list(c("Baret1991", "Normalised Ratio Vegetation Index" ), + function(red, nir) {(red/nir - 1)/(red/nir + 1)}), + REIP = list(c("GuyotAndBarnet1988", "Red Edge Inflection Point"), + function(red, redEdge1, redEdge2, redEdge3) {0.705+0.35*((red+redEdge3)/(2-redEdge1))/(redEdge2-redEdge1)}), + RVI = list(c("", "Ratio Vegetation Index"), + function(red, nir) {red/nir}), + SATVI = list(c("Marsett2006", "Soil Adjusted Total Vegetation Index"), + function(red, swir2, swir3) {(swir2 - red) / (swir2 + red + L) * (1 + L) - (swir3 / 2)}), + SAVI = list(c("Huete1988", "Soil Adjusted Vegetation Index"), + function(red, nir) {(nir - red) * (1+L) / (nir + red + L)}), + SLAVI = list(c("Lymburger2000", "Specific Leaf Area Vegetation Index"), + function(red, nir, swir2) {nir / (red + swir2)}), + SR = list(c("Birth1968", "Simple Ratio Vegetation Index"), + function(red, nir) {nir / red}), + TTVI = list(c("Thiam1997", "Thiam's Transformed Vegetation Index"), + function(red, nir) {sqrt(abs((nir-red)/(nir+red) + 0.5))}), + TVI = list(c("Deering1975", "Transformed Vegetation Index"), + function(red, nir) {sqrt((nir-red)/(nir+red)+0.5)}), + WDVI = list(c("Richardson1977", "Weighted Difference Vegetation Index"), + function(red, nir) {nir - s * red}) + ) + if(is.null(getOption("RStoolbox.verbose"))) options(RStoolbox.verbose = FALSE) + if(is.null(getOption("RStoolbox.idxdb"))) options(RStoolbox.idxdb = .IDXDB) } #' Init verbosity within functions @@ -270,6 +347,18 @@ options(RStoolbox.verbose = verbose) } + +#' Init spectralIndices within the index DB +#' +#' will restore global options after function has been called +#' @param spectralIndices List +#' @keywords internal +#' @noRd +.initIDXdb <- function(idxdb){ + idxbg <- force(getOption("RStoolbox.idxdb")) + do.call("on.exit", list(substitute(options(RStoolbox.idxdb = idxbg))), envir=parent.frame()) + options(RStoolbox.idxdb = idxdb) +} #' Clean up on package unload diff --git a/R/mesma.R b/R/mesma.R index 35e6bbc..a269399 100644 --- a/R/mesma.R +++ b/R/mesma.R @@ -10,7 +10,7 @@ #' } #' @param iterate Integer. Set maximum iteration per pixel. Processing time could increase the more iterations are made possible. Default is 400. #' @param tolerance Numeric. Tolerance limit representing a nearly zero minimal number. Default is 1e-8. -#' @param n_models Logical. Only applies if \code{em} contains column \code{class}. Defines how many endmember combinations should be picked. Maximum is the minimum number of endmembers of a class. Defaults to 5. +#' @param n_models Logical. Only applies if \code{em} contains column \code{class}. Defines how many endmember combinations should be picked. Maximum is the minimum number of endmembers of a class. Defaults to the amount of rows of em. #' @param sum_to_one Logical. Defines whether a sum-to-one constraint should be applied so that probabilities of endmember classes sum to one (a constraint not covered by NNLS) to be interpretable as fractions. Defaults to \code{TRUE}. To get actual NNLS results, change to \code{FALSE}. #' @param verbose Logical. Prints progress messages during execution. #' @param ... further arguments passed to \link[terra]{writeRaster}. @@ -72,7 +72,7 @@ #' @importFrom terra which.min rast selectRange nlyr #' @export #' -mesma <- function(img, em, method = "NNLS", iterate = 400, tolerance = 0.00000001, n_models = 5, sum_to_one = TRUE, ..., verbose){ +mesma <- function(img, em, method = "NNLS", iterate = 400, tolerance = 0.00000001, n_models = NULL, sum_to_one = TRUE, ..., verbose){ img <- .toTerra(img) ## messages if(!missing("verbose")) .initVerbose(verbose) @@ -89,6 +89,12 @@ mesma <- function(img, em, method = "NNLS", iterate = 400, tolerance = 0.0000000 if(anyNA(em)){ stop("'em' is not allowed to contain NA values. Spectra must be consistent.") } + + # Check for n_models type + if(!is.null(n_models) && !is.numeric(n_models)){ + stop("'n_models' needs to be an 'numeric' object.") + } + method <- toupper(method) meth_avail <- c("NNLS") #available methods @@ -115,6 +121,11 @@ mesma <- function(img, em, method = "NNLS", iterate = 400, tolerance = 0.0000000 classes <- unique(em$class) # check n_models n_em_cl <- min(sapply(classes, function(x) nrow(em[em$class == x,]), USE.NAMES = F)) + + if(is.null(n_models)){ + n_models <- n_em_cl + } + if(n_em_cl < n_models){ warning(paste0("'n_models' cannot be larger than the minimum number of provided endmembers per class. Setting 'n_models' to ", n_em_cl, ".")) n_models <- n_em_cl diff --git a/R/rsOpts.R b/R/rsOpts.R index aa93d38..b5e4875 100644 --- a/R/rsOpts.R +++ b/R/rsOpts.R @@ -3,14 +3,20 @@ #' shortcut to options(RStoolbox.*) #' #' @param verbose Logical. If \code{TRUE} many functions will print status messages about the current processing step. By default verbose mode is disabled. +#' @param idxdb List. The list conatins the formal calculation of spectral indices. Modify this list to pipe your own spectral index through the internal C++ calculation of RStoolbox. #' @export #' @return -#' No return, just a setter for the verbosiness of the RStoolbox package +#' No return, just a setter for the verbosiness and the index-database of the RStoolbox package. For latter, see the example of Rstoolbox::spectralIndices() #' @examples #' rsOpts(verbose=TRUE) #' -rsOpts <- function(verbose){ - options(RStoolbox.verbose=verbose) +rsOpts <- function(verbose=NULL, idxdb=NULL){ + if(!is.null(verbose)){ + options(RStoolbox.verbose=verbose) + } + if(!is.null(idxdb)){ + options(RStoolbox.idxdb=idxdb) + } } diff --git a/R/spectralIndices.R b/R/spectralIndices.R index 5af8970..13d4795 100644 --- a/R/spectralIndices.R +++ b/R/spectralIndices.R @@ -1,6 +1,6 @@ #' Spectral Indices #' -#' Calculate a suite of multispectral indices such as NDVI, SAVI etc. in an efficient way. +#' Calculate a suite of multispectral indices such as NDVI, SAVI etc. in an efficient way via C++. #' #' @param img SpatRaster. Typically remote sensing imagery, which is to be classified. #' @param blue Character or integer. Blue band. @@ -23,7 +23,7 @@ #' @param coefs List of coefficients (see Details). #' @param ... further arguments such as filename etc. passed to \link[terra]{writeRaster} #' @return SpatRaster -#' @template spectralIndices_table +#' @template spectralIndices_table #' @export #' @examples #' library(ggplot2) @@ -43,12 +43,24 @@ #' #' SI <- spectralIndices(lsat_ref, red = "B3_tre", nir = "B4_tre") #' plot(SI) +#' +#' ## Custom Spectral Index Calculation (beta) (supports only bands right now...) +#' # Get all indices +#' idxdb <- getOption("RStoolbox.idxdb") +#' +#' # Customize the RStoolbox index-database and overwrite the option +#' cdb <- c(idxdb, CUSTOM = list( list(c("Mueller2024", "Super custom index"), +#' function(blue, red) {blue + red}))) +#' rsOpts(idxdb = cdb) +#' +#' # Calculate the custom index, (also together with the provided ones) +#' custom_ind <- spectralIndices(lsat, blue = 1, red = 3, nir = 4, indices = c("NDVI", "CUSTOM")) #' } spectralIndices <- function(img, blue=NULL, green=NULL, red=NULL, nir = NULL, redEdge1 = NULL, redEdge2 = NULL, redEdge3 = NULL, swir1 = NULL, swir2 =NULL, swir3 = NULL, - scaleFactor = 1, skipRefCheck = FALSE, indices=NULL, index = NULL, maskLayer = NULL, maskValue = 1, + scaleFactor = 1, skipRefCheck = FALSE, indices=NULL, index = NULL, maskLayer = NULL, maskValue = 1, coefs = list(L = 0.5, G = 2.5, L_evi = 1, C1 = 6, C2 = 7.5, s = 1, swir2ccc = NULL, swir2coc = NULL), ... ) { # TODO: soil line estimator @@ -69,9 +81,10 @@ spectralIndices <- function(img, # mir2 | midwave infra-red | 4500 - 5000 nm # TIR1 | thermal infra-red | 8000 - 9500 nm # TIR2 | thermal infra-red | 10000 - 140000 nm - ## + ## + img <- .toTerra(img) - if(!is.null(maskLayer)) maskLayer <- .toTerra(maskLayer) + if(!is.null(maskLayer)) maskLayer <- .toTerra(maskLayer) if(!is.null(index)) indices <- index ## argument translation for convenience if("LSWI" %in% toupper(indices)) stop("LSWI has been deprecated. Use NDWI2 instead; it is identical.") if(!is.null(swir1)) stop(paste0("Currently there is no spectral index requiring swir1.", @@ -85,9 +98,13 @@ spectralIndices <- function(img, if(any(!implem)) warning("Non-implemented coefficients are ignored: ", paste0(names(coefs)[!implem], collapse=", "), "\nimplemented coefficients are: ", paste0(names(defaultCoefs), collapse = ", ")) coefs <- c(coefs, defaultCoefs[setdiff(names(defaultCoefs), names(coefs))]) - + + # Update bandsdb + idxdb <- getOption("RStoolbox.idxdb") + BANDSdb <- lapply(lapply(idxdb,"[[", 2), function(x) names(formals(x))) + ## Check indices - ind <- if(is.null(indices)) names(BANDSdb) else indices <- toupper(indices) + ind <- if(is.null(indices)) names(BANDSdb) else indices <- toupper(indices) if((is.null(coefs$swir2ccc) | is.null(coefs$swir2coc))) { if(!is.null(indices) & ("NDVIC" %in% ind)) warning("NDVIc can only be calculated if swir2ccc and swir2coc coefficients are provided.", call. = FALSE) coefs$swir2ccc <- 0 ## dummy, cant pass NULL to spectralIndicesCpp @@ -168,11 +185,9 @@ spectralIndices <- function(img, filename <- if("filename" %in% names(wopt)) wopt$filename else "" overwrite <- if("overwrite" %in% names(wopt)) wopt$overwrite else FALSE wopt[c("filename","overwrite")] <- NULL - - + # Perform calculations indexMagic <- terra::app(img[[bandsCalc]], fun = function(m) { - spectralIndicesCpp( x = m, indices = canCalc, @@ -190,7 +205,8 @@ spectralIndices <- function(img, maskValue = maskValue, L = coefs[["L"]], G = coefs[["G"]], Levi = coefs[["L_evi"]], C1 = coefs[["C1"]], C2 = coefs[["C2"]], s = coefs[["s"]], - swir2ccc = coefs[["swir2ccc"]], swir2cdiff = coefs[["swir2cdiff"]], sf = scaleFactor + swir2ccc = coefs[["swir2ccc"]], swir2cdiff = coefs[["swir2cdiff"]], sf = scaleFactor, + formulas = as.vector(unlist(lapply(idxdb[canCalc], FUN = extract_function_string))) )}, wopt=wopt, filename = filename, overwrite = overwrite) @@ -199,77 +215,21 @@ spectralIndices <- function(img, return(indexMagic) } +extract_function_string <- function(fn) { + function_body <- body(fn[[2]]) + function_string <- deparse(function_body) + function_string <- paste(function_string, collapse = "") + function_string <- gsub("\\s+", "", function_string) + function_string <- gsub("[{}]", "", function_string) + return(function_string) +} + ## NOT USED FOR CALCULATIONS ONLY FOR DOCUMENTATION ## SEE /src/spectraIndices.cpp for calculations #' Database of spectral indices #' @keywords internal -#' @noRd -.IDXdb <- list( - CLG = list(c("Gitelson2003", "Green-band Chlorophyll Index"), - function(redEdge3, green) {redEdge3/green - 1}), - CLRE = list(c("Gitelson2003", "Red-edge-band Chlorophyll Index"), - function(redEdge3, redEdge1) {redEdge3/redEdge1 - 1}), - CTVI = list(c("Perry1984", "Corrected Transformed Vegetation Index"), - function(red, nir) {(NDVI+.5)/sqrt(abs(NDVI+.5))} ), - DVI = list(c("Richardson1977", "Difference Vegetation Index"), - function(red, nir) {s*nir-red}), - EVI = list(c("Huete1999", "Enhanced Vegetation Index"), - function(red, nir, blue) {G * ((nir - red) / (nir + C1 * red - C2 * blue + L_evi))}), - EVI2 = list(c("Jiang 2008", "Two-band Enhanced Vegetation Index"), - function(red, nir) {G * (nir-red)/(nir + 2.4*red +1)}), - GEMI = list(c("Pinty1992", "Global Environmental Monitoring Index"), - function(red, nir) {(((nir^2 - red^2) * 2 + (nir * 1.5) + (red * 0.5) ) / (nir + red + 0.5)) * (1 - ((((nir^2 - red^2) * 2 + (nir * 1.5) + (red * 0.5) ) / (nir + red + 0.5)) * 0.25)) - ((red - 0.125) / (1 - red))}), - GNDVI = list(c("Gitelson1998", "Green Normalised Difference Vegetation Index" ), - function(green, nir) {(nir-green)/(nir+green)}), - KNDVI = list(c("Camps-Valls2021", "Kernel Normalised Difference Vegetation Index"), - function(red, nir) {tanh(((nir-red)/(nir+red)))^2}), - MCARI = list(c("Daughtery2000", "Modified Chlorophyll Absorption Ratio Index" ), - function(green, red, redEdge1) {((redEdge1-red)-(redEdge1-green))*(redEdge1/red)}), - MNDWI = list(c("Xu2006", "Modified Normalised Difference Water Index"), - function(green, swir2) {(green-swir2) / (green+swir2)}), - MSAVI = list(c("Qi1994", "Modified Soil Adjusted Vegetation Index" ), - function(red, nir) {nir + 0.5 - (0.5 * sqrt((2 * nir + 1)^2 - 8 * (nir - (2 * red))))}), - MSAVI2 = list(c("Qi1994", "Modified Soil Adjusted Vegetation Index 2" ), - function(red, nir) {(2 * (nir + 1) - sqrt((2 * nir + 1)^2 - 8 * (nir - red))) / 2}), - MTCI = list(c("DashAndCurran2004", "MERIS Terrestrial Chlorophyll Index"), - function(red, redEdge1, redEdge2) {(redEdge2-redEdge1)/(redEdge1-red)}), - NBRI = list(c("Garcia1991", "Normalised Burn Ratio Index"), - function(nir, swir3) { (nir - swir3) / (nir + swir3)}), - NDREI1 = list(c("GitelsonAndMerzlyak1994", "Normalised Difference Red Edge Index 1"), - function(redEdge2, redEdge1) {(redEdge2-redEdge1)/(redEdge2+redEdge1)}), - NDREI2 = list(c("Barnes2000", "Normalised Difference Red Edge Index 2"), - function(redEdge3, redEdge1) {(redEdge3-redEdge1)/(redEdge3+redEdge1)}), - NDVI = list(c("Rouse1974", "Normalised Difference Vegetation Index"), - function(red, nir) {(nir-red)/(nir+red)}), - NDVIC = list(c("Nemani1993", "Corrected Normalised Difference Vegetation Index" ), - function(red, nir, swir2) {(nir-red)/(nir+red)*(1-((swir2 - swir2ccc)/(swir2coc-swir2ccc)))}), - NDWI = list(c("McFeeters1996", "Normalised Difference Water Index"), - function(green, nir) {(green - nir)/(green + nir)}), - NDWI2 = list(c("Gao1996", "Normalised Difference Water Index"), - function(nir, swir2) {(nir - swir2)/(nir + swir2)}), - NRVI = list(c("Baret1991", "Normalised Ratio Vegetation Index" ), - function(red, nir) {(red/nir - 1)/(red/nir + 1)}), - REIP = list(c("GuyotAndBarnet1988", "Red Edge Inflection Point"), - function(red, redEdge1, redEdge2, redEdge3) {0.705+0.35*((red+redEdge3)/(2-redEdge1))/(redEdge2-redEdge1)}), - RVI = list(c("", "Ratio Vegetation Index"), - function(red, nir) {red/nir}), - SATVI = list(c("Marsett2006", "Soil Adjusted Total Vegetation Index"), - function(red, swir2, swir3) {(swir2 - red) / (swir2 + red + L) * (1 + L) - (swir3 / 2)}), - SAVI = list(c("Huete1988", "Soil Adjusted Vegetation Index"), - function(red, nir) {(nir - red) * (1+L) / (nir + red + L)}), - SLAVI = list(c("Lymburger2000", "Specific Leaf Area Vegetation Index"), - function(red, nir, swir2) {nir / (red + swir2)}), - SR = list(c("Birth1968", "Simple Ratio Vegetation Index"), - function(red, nir) {nir / red}), - TTVI = list(c("Thiam1997", "Thiam's Transformed Vegetation Index"), - function(red, nir) {sqrt(abs((nir-red)/(nir+red) + 0.5))}), - TVI = list(c("Deering1975", "Transformed Vegetation Index"), - function(red, nir) {sqrt((nir-red)/(nir+red)+0.5)}), - WDVI = list(c("Richardson1977", "Weighted Difference Vegetation Index"), - function(red, nir) {nir - s * red}) -) - - -BANDSdb <- lapply(lapply(.IDXdb,"[[", 2), function(x) names(formals(x))) +#' @noRd +.IDXdb <- getOption("RStoolbox.idxdb") +# BANDSdb <- lapply(lapply(.IDXdb,"[[", 2), function(x) names(formals(x))) .IDXdbFormulae <- lapply(.IDXdb,"[[",2) .IDX.REFdb <- lapply(.IDXdb,"[[",1) \ No newline at end of file diff --git a/R/superClass.R b/R/superClass.R index afb314b..fb0927c 100644 --- a/R/superClass.R +++ b/R/superClass.R @@ -23,7 +23,8 @@ #' @param verbose Logical. prints progress and statistics during execution #' @param overwrite logical. Overwrite spatial prediction raster if it already exists. #' @param ... further arguments to be passed to \code{\link[caret]{train}} -#' @details +#' @details +#' Note that superClass automatically loads the lattice and randomForest package. #' SuperClass performs the following steps: #' #' \enumerate{ diff --git a/man/RStoolbox.Rd b/man/RStoolbox.Rd index c725494..387edf8 100755 --- a/man/RStoolbox.Rd +++ b/man/RStoolbox.Rd @@ -64,6 +64,7 @@ The RStoolbox package provides a set of functions which simplify performing stan } } +\keyword{"RStoolbox"} \keyword{earth-observation} \keyword{remote-sensing} \keyword{spatial-data-analysis} diff --git a/man/mesma.Rd b/man/mesma.Rd index fc18ed0..35a2019 100644 --- a/man/mesma.Rd +++ b/man/mesma.Rd @@ -10,7 +10,7 @@ mesma( method = "NNLS", iterate = 400, tolerance = 1e-08, - n_models = 5, + n_models = NULL, sum_to_one = TRUE, ..., verbose @@ -30,7 +30,7 @@ mesma( \item{tolerance}{Numeric. Tolerance limit representing a nearly zero minimal number. Default is 1e-8.} -\item{n_models}{Logical. Only applies if \code{em} contains column \code{class}. Defines how many endmember combinations should be picked. Maximum is the minimum number of endmembers of a class. Defaults to 5.} +\item{n_models}{Logical. Only applies if \code{em} contains column \code{class}. Defines how many endmember combinations should be picked. Maximum is the minimum number of endmembers of a class. Defaults to the amount of rows of em.} \item{sum_to_one}{Logical. Defines whether a sum-to-one constraint should be applied so that probabilities of endmember classes sum to one (a constraint not covered by NNLS) to be interpretable as fractions. Defaults to \code{TRUE}. To get actual NNLS results, change to \code{FALSE}.} diff --git a/man/rsOpts.Rd b/man/rsOpts.Rd index 3f2a089..8c33b6e 100755 --- a/man/rsOpts.Rd +++ b/man/rsOpts.Rd @@ -4,13 +4,15 @@ \alias{rsOpts} \title{Set global options for RStoolbox} \usage{ -rsOpts(verbose) +rsOpts(verbose = NULL, idxdb = NULL) } \arguments{ \item{verbose}{Logical. If \code{TRUE} many functions will print status messages about the current processing step. By default verbose mode is disabled.} + +\item{idxdb}{List. The list conatins the formal calculation of spectral indices. Modify this list to pipe your own spectral index through the internal C++ calculation of RStoolbox.} } \value{ -No return, just a setter for the verbosiness of the RStoolbox package +No return, just a setter for the verbosiness and the index-database of the RStoolbox package. For latter, see the example of Rstoolbox::spectralIndices() } \description{ shortcut to options(RStoolbox.*) diff --git a/man/spectralIndices.Rd b/man/spectralIndices.Rd index 33cec27..80501ed 100755 --- a/man/spectralIndices.Rd +++ b/man/spectralIndices.Rd @@ -71,7 +71,7 @@ If there are invalid reflectances, e.g. clouds with reflectance > 1 this check w SpatRaster } \description{ -Calculate a suite of multispectral indices such as NDVI, SAVI etc. in an efficient way. +Calculate a suite of multispectral indices such as NDVI, SAVI etc. in an efficient way via C++. } \details{ \code{spectralIndices} calculates all indices in one go in C++, which is more efficient than calculating each index separately (for large rasters). @@ -124,5 +124,17 @@ lsat_ref <- radCor(lsat, mtlFile, method = "apref") SI <- spectralIndices(lsat_ref, red = "B3_tre", nir = "B4_tre") plot(SI) + +## Custom Spectral Index Calculation (beta) (supports only bands right now...) +# Get all indices +idxdb <- getOption("RStoolbox.idxdb") + +# Cutomize the RStoolbox index-database and overwrite the option +cdb <- c(idxdb, CUSTOM = list( list(c("Mueller2024", "Super custom index"), + function(blue, red) {blue + red}))) +rsOpts(idxdb = cdb) + +# Calculate the custom index, (also together with the provided ones) +custom_ind <- spectralIndices(lsat, blue = 1, red = 3, nir = 4, indices = c("NDVI", "CUSTOM")) } } diff --git a/man/superClass.Rd b/man/superClass.Rd index 4a1f07f..6e9fa39 100755 --- a/man/superClass.Rd +++ b/man/superClass.Rd @@ -87,6 +87,7 @@ A superClass object (effectively a list) containing: Supervised classification both for classification and regression mode based on vector training data (points or polygons). } \details{ +Note that superClass automatically loads the lattice and randomForest package. SuperClass performs the following steps: \enumerate{ diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index cb4993f..442f8ee 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -167,8 +167,8 @@ BEGIN_RCPP END_RCPP } // spectralIndicesCpp -NumericMatrix spectralIndicesCpp(NumericMatrix x, CharacterVector indices, const int redBand, const int blueBand, const int greenBand, const int nirBand, const int redEdge1Band, const int redEdge2Band, const int redEdge3Band, const int swir1Band, const int swir2Band, const int swir3Band, int maskLayer, const int maskValue, const double L, const double s, const double G, const double C1, const double C2, double Levi, const double swir2ccc, const double swir2cdiff, const double sf); -RcppExport SEXP _RStoolbox_spectralIndicesCpp(SEXP xSEXP, SEXP indicesSEXP, SEXP redBandSEXP, SEXP blueBandSEXP, SEXP greenBandSEXP, SEXP nirBandSEXP, SEXP redEdge1BandSEXP, SEXP redEdge2BandSEXP, SEXP redEdge3BandSEXP, SEXP swir1BandSEXP, SEXP swir2BandSEXP, SEXP swir3BandSEXP, SEXP maskLayerSEXP, SEXP maskValueSEXP, SEXP LSEXP, SEXP sSEXP, SEXP GSEXP, SEXP C1SEXP, SEXP C2SEXP, SEXP LeviSEXP, SEXP swir2cccSEXP, SEXP swir2cdiffSEXP, SEXP sfSEXP) { +NumericMatrix spectralIndicesCpp(NumericMatrix x, CharacterVector indices, const int redBand, const int blueBand, const int greenBand, const int nirBand, const int redEdge1Band, const int redEdge2Band, const int redEdge3Band, const int swir1Band, const int swir2Band, const int swir3Band, int maskLayer, const int maskValue, const double L, const double s, const double G, const double C1, const double C2, double Levi, const double swir2ccc, const double swir2cdiff, const double sf, CharacterVector formulas); +RcppExport SEXP _RStoolbox_spectralIndicesCpp(SEXP xSEXP, SEXP indicesSEXP, SEXP redBandSEXP, SEXP blueBandSEXP, SEXP greenBandSEXP, SEXP nirBandSEXP, SEXP redEdge1BandSEXP, SEXP redEdge2BandSEXP, SEXP redEdge3BandSEXP, SEXP swir1BandSEXP, SEXP swir2BandSEXP, SEXP swir3BandSEXP, SEXP maskLayerSEXP, SEXP maskValueSEXP, SEXP LSEXP, SEXP sSEXP, SEXP GSEXP, SEXP C1SEXP, SEXP C2SEXP, SEXP LeviSEXP, SEXP swir2cccSEXP, SEXP swir2cdiffSEXP, SEXP sfSEXP, SEXP formulasSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -195,7 +195,8 @@ BEGIN_RCPP Rcpp::traits::input_parameter< const double >::type swir2ccc(swir2cccSEXP); Rcpp::traits::input_parameter< const double >::type swir2cdiff(swir2cdiffSEXP); Rcpp::traits::input_parameter< const double >::type sf(sfSEXP); - rcpp_result_gen = Rcpp::wrap(spectralIndicesCpp(x, indices, redBand, blueBand, greenBand, nirBand, redEdge1Band, redEdge2Band, redEdge3Band, swir1Band, swir2Band, swir3Band, maskLayer, maskValue, L, s, G, C1, C2, Levi, swir2ccc, swir2cdiff, sf)); + Rcpp::traits::input_parameter< CharacterVector >::type formulas(formulasSEXP); + rcpp_result_gen = Rcpp::wrap(spectralIndicesCpp(x, indices, redBand, blueBand, greenBand, nirBand, redEdge1Band, redEdge2Band, redEdge3Band, swir1Band, swir2Band, swir3Band, maskLayer, maskValue, L, s, G, C1, C2, Levi, swir2ccc, swir2cdiff, sf, formulas)); return rcpp_result_gen; END_RCPP } diff --git a/src/init.c b/src/init.c index 0b015ef..5afa3c7 100644 --- a/src/init.c +++ b/src/init.c @@ -19,7 +19,7 @@ extern SEXP _RStoolbox_predKmeansCpp(SEXP, SEXP, SEXP); extern SEXP _RStoolbox_pwSimilarityCpp(SEXP, SEXP, SEXP); extern SEXP _RStoolbox_rescaleImageCpp(SEXP, SEXP, SEXP, SEXP); extern SEXP _RStoolbox_specSimC(SEXP, SEXP); -extern SEXP _RStoolbox_spectralIndicesCpp(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +extern SEXP _RStoolbox_spectralIndicesCpp(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); extern SEXP _RStoolbox_availableRAMCpp(SEXP); static const R_CallMethodDef CallEntries[] = { @@ -34,7 +34,7 @@ static const R_CallMethodDef CallEntries[] = { {"_RStoolbox_pwSimilarityCpp", (DL_FUNC) &_RStoolbox_pwSimilarityCpp, 3}, {"_RStoolbox_rescaleImageCpp", (DL_FUNC) &_RStoolbox_rescaleImageCpp, 4}, {"_RStoolbox_specSimC", (DL_FUNC) &_RStoolbox_specSimC, 2}, - {"_RStoolbox_spectralIndicesCpp", (DL_FUNC) &_RStoolbox_spectralIndicesCpp, 23}, + {"_RStoolbox_spectralIndicesCpp", (DL_FUNC) &_RStoolbox_spectralIndicesCpp, 24}, {"_RStoolbox_availableRAMCpp", (DL_FUNC) &_RStoolbox_availableRAMCpp, 1}, {NULL, NULL, 0} }; diff --git a/src/spectralIndices.cpp b/src/spectralIndices.cpp index c9f6674..5e89d00 100644 --- a/src/spectralIndices.cpp +++ b/src/spectralIndices.cpp @@ -1,6 +1,10 @@ #include +#include +#include +#include "tinyexpr.h" using namespace Rcpp; +using namespace std; //[[Rcpp::export]] NumericMatrix spectralIndicesCpp(NumericMatrix x, CharacterVector indices, @@ -9,7 +13,8 @@ NumericMatrix spectralIndicesCpp(NumericMatrix x, CharacterVector indices, const int swir1Band, const int swir2Band, const int swir3Band, int maskLayer, const int maskValue, const double L, const double s, const double G, const double C1, - const double C2, double Levi, const double swir2ccc, const double swir2cdiff, const double sf) { + const double C2, double Levi, const double swir2ccc, const double swir2cdiff, const double sf, + CharacterVector formulas) { const int nind = indices.size(); const int nsamp = x.nrow(); @@ -17,7 +22,7 @@ NumericMatrix spectralIndicesCpp(NumericMatrix x, CharacterVector indices, NumericMatrix out(nsamp, nind); NumericVector blue, green, red, redEdge1, redEdge2, redEdge3, nir, swir1, swir2, swir3; - + // Apply mask layer if(!IntegerVector::is_na(maskLayer)){ maskLayer-=1 ; @@ -235,6 +240,83 @@ NumericMatrix spectralIndicesCpp(NumericMatrix x, CharacterVector indices, out(_,j) = nir - s * red; out(_,j) = ifelse(is_na(out(_,j)), NA_REAL, out(_,j)); } + else { + std::string formula = Rcpp::as(formulas[j]); + //std::string formula = "red * red"; + + size_t vec_size = 0; + bool red_present = formula.find("red") != std::string::npos; + bool blue_present = formula.find("blue") != std::string::npos; + bool green_present = formula.find("green") != std::string::npos; + bool redEdge1_present = formula.find("redEdge1") != std::string::npos; + bool redEdge2_present = formula.find("redEdge2") != std::string::npos; + bool redEdge3_present = formula.find("redEdge3") != std::string::npos; + bool nir_present = formula.find("nir") != std::string::npos; + bool swir1_present = formula.find("swir1") != std::string::npos; + bool swir2_present = formula.find("swir2") != std::string::npos; + bool swir3_present = formula.find("swir3") != std::string::npos; + + if (red_present) vec_size = red.size(); + else if (blue_present) vec_size = blue.size(); + else if (green_present) vec_size = green.size(); + else if (redEdge1_present) vec_size = redEdge1.size(); + else if (redEdge2_present) vec_size = redEdge2.size(); + else if (redEdge3_present) vec_size = redEdge3.size(); + else if (nir_present) vec_size = nir.size(); + else if (swir1_present) vec_size = swir1.size(); + else if (swir2_present) vec_size = swir2.size(); + else if (swir3_present) vec_size = swir3.size(); + else { + Rcpp::stop("No valid variable found in the formula."); + } + + std::vector result(vec_size); + + std::vector vars; + + std::vector valid_var_names = { + "blue", "green", "red", "redEdge1", "redEdge2", "redEdge3", "nir", "swri1", "swir2", "swir3" + }; + std::vector var_stores(valid_var_names.size()); + + for (size_t s = 0; s < valid_var_names.size(); ++s) { + if(formula.find(valid_var_names[s]) != std::string::npos){ + vars.push_back({valid_var_names[s].c_str(), &var_stores[s]}); + } + } + + int err; + te_expr* expr = te_compile(formula.c_str(), vars.data(), vars.size(), &err); + + if (expr) { + try { + for (size_t i = 0; i < vec_size; ++i) { + if (red_present) var_stores[std::distance(valid_var_names.begin(), std::find(valid_var_names.begin(), valid_var_names.end(), "red"))] = red[i]; + if (blue_present) var_stores[std::distance(valid_var_names.begin(), std::find(valid_var_names.begin(), valid_var_names.end(), "blue"))] = blue[i]; + if (green_present) var_stores[std::distance(valid_var_names.begin(), std::find(valid_var_names.begin(), valid_var_names.end(), "green"))] = green[i]; + if (redEdge1_present) var_stores[std::distance(valid_var_names.begin(), std::find(valid_var_names.begin(), valid_var_names.end(), "redEdge1"))] = redEdge1[i]; + if (redEdge2_present) var_stores[std::distance(valid_var_names.begin(), std::find(valid_var_names.begin(), valid_var_names.end(), "redEdge2"))] = redEdge2[i]; + if (redEdge3_present) var_stores[std::distance(valid_var_names.begin(), std::find(valid_var_names.begin(), valid_var_names.end(), "redEdge3"))] = redEdge3[i]; + if (nir_present) var_stores[std::distance(valid_var_names.begin(), std::find(valid_var_names.begin(), valid_var_names.end(), "nir"))] = nir[i]; + if (swir1_present) var_stores[std::distance(valid_var_names.begin(), std::find(valid_var_names.begin(), valid_var_names.end(), "swir1"))] = swir1[i]; + if (swir2_present) var_stores[std::distance(valid_var_names.begin(), std::find(valid_var_names.begin(), valid_var_names.end(), "swir2"))] = swir2[i]; + if (swir3_present) var_stores[std::distance(valid_var_names.begin(), std::find(valid_var_names.begin(), valid_var_names.end(), "swir3"))] = swir3[i]; + + result[i] = te_eval(expr); + } + } catch (...) { + Rcpp::stop("Error occurred during evaluation. Ensure all required variables are provided."); + } + + te_free(expr); + } else { + Rcpp::stop("Failed to parse the formula. Error code: " + std::to_string(err)); + } + + Rcpp::NumericVector r_result(result.begin(), result.end()); + out(_,j) = r_result; + + } } diff --git a/src/tinyexpr.c b/src/tinyexpr.c new file mode 100644 index 0000000..3e7ba4f --- /dev/null +++ b/src/tinyexpr.c @@ -0,0 +1,705 @@ +// SPDX-License-Identifier: Zlib +/* + * TINYEXPR - Tiny recursive descent parser and evaluation engine in C + * + * Copyright (c) 2015-2020 Lewis Van Winkle + * + * http://CodePlea.com + * + * This software is provided 'as-is', without any express or implied + * warranty. In no event will the authors be held liable for any damages + * arising from the use of this software. + * + * Permission is granted to anyone to use this software for any purpose, + * including commercial applications, and to alter it and redistribute it + * freely, subject to the following restrictions: + * + * 1. The origin of this software must not be misrepresented; you must not + * claim that you wrote the original software. If you use this software + * in a product, an acknowledgement in the product documentation would be + * appreciated but is not required. + * 2. Altered source versions must be plainly marked as such, and must not be + * misrepresented as being the original software. + * 3. This notice may not be removed or altered from any source distribution. + */ + +/* COMPILE TIME OPTIONS */ + +/* Exponentiation associativity: +For a^b^c = (a^b)^c and -a^b = (-a)^b do nothing. +For a^b^c = a^(b^c) and -a^b = -(a^b) uncomment the next line.*/ +/* #define TE_POW_FROM_RIGHT */ + +/* Logarithms +For log = base 10 log do nothing +For log = natural log uncomment the next line. */ +/* #define TE_NAT_LOG */ + +#include "tinyexpr.h" +#include +#include +#include +#include +#include +#include + +#ifndef NAN +#define NAN (0.0/0.0) +#endif + +#ifndef INFINITY +#define INFINITY (1.0/0.0) +#endif + + +typedef double (*te_fun2)(double, double); + +enum { + TOK_NULL = TE_CLOSURE7+1, TOK_ERROR, TOK_END, TOK_SEP, + TOK_OPEN, TOK_CLOSE, TOK_NUMBER, TOK_VARIABLE, TOK_INFIX +}; + + +enum {TE_CONSTANT = 1}; + + +typedef struct state { + const char *start; + const char *next; + int type; + union {double value; const double *bound; const void *function;}; + void *context; + + const te_variable *lookup; + int lookup_len; +} state; + + +#define TYPE_MASK(TYPE) ((TYPE)&0x0000001F) + +#define IS_PURE(TYPE) (((TYPE) & TE_FLAG_PURE) != 0) +#define IS_FUNCTION(TYPE) (((TYPE) & TE_FUNCTION0) != 0) +#define IS_CLOSURE(TYPE) (((TYPE) & TE_CLOSURE0) != 0) +#define ARITY(TYPE) ( ((TYPE) & (TE_FUNCTION0 | TE_CLOSURE0)) ? ((TYPE) & 0x00000007) : 0 ) +#define NEW_EXPR(type, ...) new_expr((type), (const te_expr*[]){__VA_ARGS__}) +#define CHECK_NULL(ptr, ...) if ((ptr) == NULL) { __VA_ARGS__; return NULL; } + +static te_expr *new_expr(const int type, const te_expr *parameters[]) { + const int arity = ARITY(type); + const int psize = sizeof(void*) * arity; + const int size = sizeof(te_expr) + psize + (IS_CLOSURE(type) ? sizeof(void*) : 0); + te_expr *ret = malloc(size); + + CHECK_NULL(ret); + + memset(ret, 0, size); + if (arity && parameters) { + memcpy(ret->parameters, parameters, psize); + } + ret->type = type; + ret->bound = 0; + return ret; +} + + +void te_free_parameters(te_expr *n) { + if (!n) return; + switch (TYPE_MASK(n->type)) { + case TE_FUNCTION7: case TE_CLOSURE7: te_free(n->parameters[6]); /* Falls through. */ + case TE_FUNCTION6: case TE_CLOSURE6: te_free(n->parameters[5]); /* Falls through. */ + case TE_FUNCTION5: case TE_CLOSURE5: te_free(n->parameters[4]); /* Falls through. */ + case TE_FUNCTION4: case TE_CLOSURE4: te_free(n->parameters[3]); /* Falls through. */ + case TE_FUNCTION3: case TE_CLOSURE3: te_free(n->parameters[2]); /* Falls through. */ + case TE_FUNCTION2: case TE_CLOSURE2: te_free(n->parameters[1]); /* Falls through. */ + case TE_FUNCTION1: case TE_CLOSURE1: te_free(n->parameters[0]); + } +} + + +void te_free(te_expr *n) { + if (!n) return; + te_free_parameters(n); + free(n); +} + + +static double pi(void) {return 3.14159265358979323846;} +static double e(void) {return 2.71828182845904523536;} +static double fac(double a) {/* simplest version of fac */ + if (a < 0.0) + return NAN; + if (a > UINT_MAX) + return INFINITY; + unsigned int ua = (unsigned int)(a); + unsigned long int result = 1, i; + for (i = 1; i <= ua; i++) { + if (i > ULONG_MAX / result) + return INFINITY; + result *= i; + } + return (double)result; +} +static double ncr(double n, double r) { + if (n < 0.0 || r < 0.0 || n < r) return NAN; + if (n > UINT_MAX || r > UINT_MAX) return INFINITY; + unsigned long int un = (unsigned int)(n), ur = (unsigned int)(r), i; + unsigned long int result = 1; + if (ur > un / 2) ur = un - ur; + for (i = 1; i <= ur; i++) { + if (result > ULONG_MAX / (un - ur + i)) + return INFINITY; + result *= un - ur + i; + result /= i; + } + return result; +} +static double npr(double n, double r) {return ncr(n, r) * fac(r);} + +#ifdef _MSC_VER +#pragma function (ceil) +#pragma function (floor) +#endif + +static const te_variable functions[] = { + /* must be in alphabetical order */ + {"abs", fabs, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"acos", acos, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"asin", asin, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"atan", atan, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"atan2", atan2, TE_FUNCTION2 | TE_FLAG_PURE, 0}, + {"ceil", ceil, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"cos", cos, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"cosh", cosh, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"e", e, TE_FUNCTION0 | TE_FLAG_PURE, 0}, + {"exp", exp, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"fac", fac, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"floor", floor, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"ln", log, TE_FUNCTION1 | TE_FLAG_PURE, 0}, +#ifdef TE_NAT_LOG + {"log", log, TE_FUNCTION1 | TE_FLAG_PURE, 0}, +#else + {"log", log10, TE_FUNCTION1 | TE_FLAG_PURE, 0}, +#endif + {"log10", log10, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"ncr", ncr, TE_FUNCTION2 | TE_FLAG_PURE, 0}, + {"npr", npr, TE_FUNCTION2 | TE_FLAG_PURE, 0}, + {"pi", pi, TE_FUNCTION0 | TE_FLAG_PURE, 0}, + {"pow", pow, TE_FUNCTION2 | TE_FLAG_PURE, 0}, + {"sin", sin, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"sinh", sinh, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"sqrt", sqrt, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"tan", tan, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"tanh", tanh, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {0, 0, 0, 0} +}; + +static const te_variable *find_builtin(const char *name, int len) { + int imin = 0; + int imax = sizeof(functions) / sizeof(te_variable) - 2; + + /*Binary search.*/ + while (imax >= imin) { + const int i = (imin + ((imax-imin)/2)); + int c = strncmp(name, functions[i].name, len); + if (!c) c = '\0' - functions[i].name[len]; + if (c == 0) { + return functions + i; + } else if (c > 0) { + imin = i + 1; + } else { + imax = i - 1; + } + } + + return 0; +} + +static const te_variable *find_lookup(const state *s, const char *name, int len) { + int iters; + const te_variable *var; + if (!s->lookup) return 0; + + for (var = s->lookup, iters = s->lookup_len; iters; ++var, --iters) { + if (strncmp(name, var->name, len) == 0 && var->name[len] == '\0') { + return var; + } + } + return 0; +} + + + +static double add(double a, double b) {return a + b;} +static double sub(double a, double b) {return a - b;} +static double mul(double a, double b) {return a * b;} +static double divide(double a, double b) {return a / b;} +static double negate(double a) {return -a;} +static double comma(double a, double b) {(void)a; return b;} + + +void next_token(state *s) { + s->type = TOK_NULL; + + do { + + if (!*s->next){ + s->type = TOK_END; + return; + } + + /* Try reading a number. */ + if ((s->next[0] >= '0' && s->next[0] <= '9') || s->next[0] == '.') { + s->value = strtod(s->next, (char**)&s->next); + s->type = TOK_NUMBER; + } else { + /* Look for a variable or builtin function call. */ + if (isalpha(s->next[0])) { + const char *start; + start = s->next; + while (isalpha(s->next[0]) || isdigit(s->next[0]) || (s->next[0] == '_')) s->next++; + + const te_variable *var = find_lookup(s, start, s->next - start); + if (!var) var = find_builtin(start, s->next - start); + + if (!var) { + s->type = TOK_ERROR; + } else { + switch(TYPE_MASK(var->type)) + { + case TE_VARIABLE: + s->type = TOK_VARIABLE; + s->bound = var->address; + break; + + case TE_CLOSURE0: case TE_CLOSURE1: case TE_CLOSURE2: case TE_CLOSURE3: /* Falls through. */ + case TE_CLOSURE4: case TE_CLOSURE5: case TE_CLOSURE6: case TE_CLOSURE7: /* Falls through. */ + s->context = var->context; /* Falls through. */ + + case TE_FUNCTION0: case TE_FUNCTION1: case TE_FUNCTION2: case TE_FUNCTION3: /* Falls through. */ + case TE_FUNCTION4: case TE_FUNCTION5: case TE_FUNCTION6: case TE_FUNCTION7: /* Falls through. */ + s->type = var->type; + s->function = var->address; + break; + } + } + + } else { + /* Look for an operator or special character. */ + switch (s->next++[0]) { + case '+': s->type = TOK_INFIX; s->function = add; break; + case '-': s->type = TOK_INFIX; s->function = sub; break; + case '*': s->type = TOK_INFIX; s->function = mul; break; + case '/': s->type = TOK_INFIX; s->function = divide; break; + case '^': s->type = TOK_INFIX; s->function = pow; break; + case '%': s->type = TOK_INFIX; s->function = fmod; break; + case '(': s->type = TOK_OPEN; break; + case ')': s->type = TOK_CLOSE; break; + case ',': s->type = TOK_SEP; break; + case ' ': case '\t': case '\n': case '\r': break; + default: s->type = TOK_ERROR; break; + } + } + } + } while (s->type == TOK_NULL); +} + + +static te_expr *list(state *s); +static te_expr *expr(state *s); +static te_expr *power(state *s); + +static te_expr *base(state *s) { + /* = | | {"(" ")"} | | "(" {"," } ")" | "(" ")" */ + te_expr *ret; + int arity; + + switch (TYPE_MASK(s->type)) { + case TOK_NUMBER: + ret = new_expr(TE_CONSTANT, 0); + CHECK_NULL(ret); + + ret->value = s->value; + next_token(s); + break; + + case TOK_VARIABLE: + ret = new_expr(TE_VARIABLE, 0); + CHECK_NULL(ret); + + ret->bound = s->bound; + next_token(s); + break; + + case TE_FUNCTION0: + case TE_CLOSURE0: + ret = new_expr(s->type, 0); + CHECK_NULL(ret); + + ret->function = s->function; + if (IS_CLOSURE(s->type)) ret->parameters[0] = s->context; + next_token(s); + if (s->type == TOK_OPEN) { + next_token(s); + if (s->type != TOK_CLOSE) { + s->type = TOK_ERROR; + } else { + next_token(s); + } + } + break; + + case TE_FUNCTION1: + case TE_CLOSURE1: + ret = new_expr(s->type, 0); + CHECK_NULL(ret); + + ret->function = s->function; + if (IS_CLOSURE(s->type)) ret->parameters[1] = s->context; + next_token(s); + ret->parameters[0] = power(s); + CHECK_NULL(ret->parameters[0], te_free(ret)); + break; + + case TE_FUNCTION2: case TE_FUNCTION3: case TE_FUNCTION4: + case TE_FUNCTION5: case TE_FUNCTION6: case TE_FUNCTION7: + case TE_CLOSURE2: case TE_CLOSURE3: case TE_CLOSURE4: + case TE_CLOSURE5: case TE_CLOSURE6: case TE_CLOSURE7: + arity = ARITY(s->type); + + ret = new_expr(s->type, 0); + CHECK_NULL(ret); + + ret->function = s->function; + if (IS_CLOSURE(s->type)) ret->parameters[arity] = s->context; + next_token(s); + + if (s->type != TOK_OPEN) { + s->type = TOK_ERROR; + } else { + int i; + for(i = 0; i < arity; i++) { + next_token(s); + ret->parameters[i] = expr(s); + CHECK_NULL(ret->parameters[i], te_free(ret)); + + if(s->type != TOK_SEP) { + break; + } + } + if(s->type != TOK_CLOSE || i != arity - 1) { + s->type = TOK_ERROR; + } else { + next_token(s); + } + } + + break; + + case TOK_OPEN: + next_token(s); + ret = list(s); + CHECK_NULL(ret); + + if (s->type != TOK_CLOSE) { + s->type = TOK_ERROR; + } else { + next_token(s); + } + break; + + default: + ret = new_expr(0, 0); + CHECK_NULL(ret); + + s->type = TOK_ERROR; + ret->value = NAN; + break; + } + + return ret; +} + + +static te_expr *power(state *s) { + /* = {("-" | "+")} */ + int sign = 1; + while (s->type == TOK_INFIX && (s->function == add || s->function == sub)) { + if (s->function == sub) sign = -sign; + next_token(s); + } + + te_expr *ret; + + if (sign == 1) { + ret = base(s); + } else { + te_expr *b = base(s); + CHECK_NULL(b); + + ret = NEW_EXPR(TE_FUNCTION1 | TE_FLAG_PURE, b); + CHECK_NULL(ret, te_free(b)); + + ret->function = negate; + } + + return ret; +} + +#ifdef TE_POW_FROM_RIGHT +static te_expr *factor(state *s) { + /* = {"^" } */ + te_expr *ret = power(s); + CHECK_NULL(ret); + + int neg = 0; + + if (ret->type == (TE_FUNCTION1 | TE_FLAG_PURE) && ret->function == negate) { + te_expr *se = ret->parameters[0]; + free(ret); + ret = se; + neg = 1; + } + + te_expr *insertion = 0; + + while (s->type == TOK_INFIX && (s->function == pow)) { + te_fun2 t = s->function; + next_token(s); + + if (insertion) { + /* Make exponentiation go right-to-left. */ + te_expr *p = power(s); + CHECK_NULL(p, te_free(ret)); + + te_expr *insert = NEW_EXPR(TE_FUNCTION2 | TE_FLAG_PURE, insertion->parameters[1], p); + CHECK_NULL(insert, te_free(p), te_free(ret)); + + insert->function = t; + insertion->parameters[1] = insert; + insertion = insert; + } else { + te_expr *p = power(s); + CHECK_NULL(p, te_free(ret)); + + te_expr *prev = ret; + ret = NEW_EXPR(TE_FUNCTION2 | TE_FLAG_PURE, ret, p); + CHECK_NULL(ret, te_free(p), te_free(prev)); + + ret->function = t; + insertion = ret; + } + } + + if (neg) { + te_expr *prev = ret; + ret = NEW_EXPR(TE_FUNCTION1 | TE_FLAG_PURE, ret); + CHECK_NULL(ret, te_free(prev)); + + ret->function = negate; + } + + return ret; +} +#else +static te_expr *factor(state *s) { + /* = {"^" } */ + te_expr *ret = power(s); + CHECK_NULL(ret); + + while (s->type == TOK_INFIX && (s->function == pow)) { + te_fun2 t = s->function; + next_token(s); + te_expr *p = power(s); + CHECK_NULL(p, te_free(ret)); + + te_expr *prev = ret; + ret = NEW_EXPR(TE_FUNCTION2 | TE_FLAG_PURE, ret, p); + CHECK_NULL(ret, te_free(p), te_free(prev)); + + ret->function = t; + } + + return ret; +} +#endif + + + +static te_expr *term(state *s) { + /* = {("*" | "/" | "%") } */ + te_expr *ret = factor(s); + CHECK_NULL(ret); + + while (s->type == TOK_INFIX && (s->function == mul || s->function == divide || s->function == fmod)) { + te_fun2 t = s->function; + next_token(s); + te_expr *f = factor(s); + CHECK_NULL(f, te_free(ret)); + + te_expr *prev = ret; + ret = NEW_EXPR(TE_FUNCTION2 | TE_FLAG_PURE, ret, f); + CHECK_NULL(ret, te_free(f), te_free(prev)); + + ret->function = t; + } + + return ret; +} + + +static te_expr *expr(state *s) { + /* = {("+" | "-") } */ + te_expr *ret = term(s); + CHECK_NULL(ret); + + while (s->type == TOK_INFIX && (s->function == add || s->function == sub)) { + te_fun2 t = s->function; + next_token(s); + te_expr *te = term(s); + CHECK_NULL(te, te_free(ret)); + + te_expr *prev = ret; + ret = NEW_EXPR(TE_FUNCTION2 | TE_FLAG_PURE, ret, te); + CHECK_NULL(ret, te_free(te), te_free(prev)); + + ret->function = t; + } + + return ret; +} + + +static te_expr *list(state *s) { + /* = {"," } */ + te_expr *ret = expr(s); + CHECK_NULL(ret); + + while (s->type == TOK_SEP) { + next_token(s); + te_expr *e = expr(s); + CHECK_NULL(e, te_free(ret)); + + te_expr *prev = ret; + ret = NEW_EXPR(TE_FUNCTION2 | TE_FLAG_PURE, ret, e); + CHECK_NULL(ret, te_free(e), te_free(prev)); + + ret->function = comma; + } + + return ret; +} + + +#define TE_FUN(...) ((double(*)(__VA_ARGS__))n->function) +#define M(e) te_eval(n->parameters[e]) + + +double te_eval(const te_expr *n) { + if (!n) return NAN; + + switch(TYPE_MASK(n->type)) { + case TE_CONSTANT: return n->value; + case TE_VARIABLE: return *n->bound; + + case TE_FUNCTION0: case TE_FUNCTION1: case TE_FUNCTION2: case TE_FUNCTION3: + case TE_FUNCTION4: case TE_FUNCTION5: case TE_FUNCTION6: case TE_FUNCTION7: + switch(ARITY(n->type)) { + case 0: return TE_FUN(void)(); + case 1: return TE_FUN(double)(M(0)); + case 2: return TE_FUN(double, double)(M(0), M(1)); + case 3: return TE_FUN(double, double, double)(M(0), M(1), M(2)); + case 4: return TE_FUN(double, double, double, double)(M(0), M(1), M(2), M(3)); + case 5: return TE_FUN(double, double, double, double, double)(M(0), M(1), M(2), M(3), M(4)); + case 6: return TE_FUN(double, double, double, double, double, double)(M(0), M(1), M(2), M(3), M(4), M(5)); + case 7: return TE_FUN(double, double, double, double, double, double, double)(M(0), M(1), M(2), M(3), M(4), M(5), M(6)); + default: return NAN; + } + + case TE_CLOSURE0: case TE_CLOSURE1: case TE_CLOSURE2: case TE_CLOSURE3: + case TE_CLOSURE4: case TE_CLOSURE5: case TE_CLOSURE6: case TE_CLOSURE7: + switch(ARITY(n->type)) { + case 0: return TE_FUN(void*)(n->parameters[0]); + case 1: return TE_FUN(void*, double)(n->parameters[1], M(0)); + case 2: return TE_FUN(void*, double, double)(n->parameters[2], M(0), M(1)); + case 3: return TE_FUN(void*, double, double, double)(n->parameters[3], M(0), M(1), M(2)); + case 4: return TE_FUN(void*, double, double, double, double)(n->parameters[4], M(0), M(1), M(2), M(3)); + case 5: return TE_FUN(void*, double, double, double, double, double)(n->parameters[5], M(0), M(1), M(2), M(3), M(4)); + case 6: return TE_FUN(void*, double, double, double, double, double, double)(n->parameters[6], M(0), M(1), M(2), M(3), M(4), M(5)); + case 7: return TE_FUN(void*, double, double, double, double, double, double, double)(n->parameters[7], M(0), M(1), M(2), M(3), M(4), M(5), M(6)); + default: return NAN; + } + + default: return NAN; + } + +} + +#undef TE_FUN +#undef M + +static void optimize(te_expr *n) { + /* Evaluates as much as possible. */ + if (n->type == TE_CONSTANT) return; + if (n->type == TE_VARIABLE) return; + + /* Only optimize out functions flagged as pure. */ + if (IS_PURE(n->type)) { + const int arity = ARITY(n->type); + int known = 1; + int i; + for (i = 0; i < arity; ++i) { + optimize(n->parameters[i]); + if (((te_expr*)(n->parameters[i]))->type != TE_CONSTANT) { + known = 0; + } + } + if (known) { + const double value = te_eval(n); + te_free_parameters(n); + n->type = TE_CONSTANT; + n->value = value; + } + } +} + + +te_expr *te_compile(const char *expression, const te_variable *variables, int var_count, int *error) { + state s; + s.start = s.next = expression; + s.lookup = variables; + s.lookup_len = var_count; + + next_token(&s); + te_expr *root = list(&s); + if (root == NULL) { + if (error) *error = -1; + return NULL; + } + + if (s.type != TOK_END) { + te_free(root); + if (error) { + *error = (s.next - s.start); + if (*error == 0) *error = 1; + } + return 0; + } else { + optimize(root); + if (error) *error = 0; + return root; + } +} + + +double te_interp(const char *expression, int *error) { + te_expr *n = te_compile(expression, 0, 0, error); + + double ret; + if (n) { + ret = te_eval(n); + te_free(n); + } else { + ret = NAN; + } + return ret; +} diff --git a/src/tinyexpr.h b/src/tinyexpr.h new file mode 100644 index 0000000..a826c8b --- /dev/null +++ b/src/tinyexpr.h @@ -0,0 +1,87 @@ +// SPDX-License-Identifier: Zlib +/* + * TINYEXPR - Tiny recursive descent parser and evaluation engine in C + * + * Copyright (c) 2015-2020 Lewis Van Winkle + * + * http://CodePlea.com + * + * This software is provided 'as-is', without any express or implied + * warranty. In no event will the authors be held liable for any damages + * arising from the use of this software. + * + * Permission is granted to anyone to use this software for any purpose, + * including commercial applications, and to alter it and redistribute it + * freely, subject to the following restrictions: + * + * 1. The origin of this software must not be misrepresented; you must not + * claim that you wrote the original software. If you use this software + * in a product, an acknowledgement in the product documentation would be + * appreciated but is not required. + * 2. Altered source versions must be plainly marked as such, and must not be + * misrepresented as being the original software. + * 3. This notice may not be removed or altered from any source distribution. + */ + +#ifndef TINYEXPR_H +#define TINYEXPR_H + + +#ifdef __cplusplus +extern "C" { +#endif + + + +typedef struct te_expr { + int type; + union {double value; const double *bound; const void *function;}; +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wpedantic" + void *parameters[]; +#pragma GCC diagnostic pop +} te_expr; + + +enum { + TE_VARIABLE = 0, + + TE_FUNCTION0 = 8, TE_FUNCTION1, TE_FUNCTION2, TE_FUNCTION3, + TE_FUNCTION4, TE_FUNCTION5, TE_FUNCTION6, TE_FUNCTION7, + + TE_CLOSURE0 = 16, TE_CLOSURE1, TE_CLOSURE2, TE_CLOSURE3, + TE_CLOSURE4, TE_CLOSURE5, TE_CLOSURE6, TE_CLOSURE7, + + TE_FLAG_PURE = 32 +}; + +typedef struct te_variable { + const char *name; + const void *address; + int type; + void *context; +} te_variable; + + + +/* Parses the input expression, evaluates it, and frees it. */ +/* Returns NaN on error. */ +double te_interp(const char *expression, int *error); + +/* Parses the input expression and binds variables. */ +/* Returns NULL on error. */ +te_expr *te_compile(const char *expression, const te_variable *variables, int var_count, int *error); + +/* Evaluates the expression. */ +double te_eval(const te_expr *n); + +/* Frees the expression. */ +/* This is safe to call on NULL pointers. */ +void te_free(te_expr *n); + + +#ifdef __cplusplus +} +#endif + +#endif /*TINYEXPR_H*/ diff --git a/tests/testthat/test-spectralIndices.R b/tests/testthat/test-spectralIndices.R index af20d17..8357797 100644 --- a/tests/testthat/test-spectralIndices.R +++ b/tests/testthat/test-spectralIndices.R @@ -69,7 +69,6 @@ test_that("returned classes", { test_that("excercise all indices", { expect_is(sp <- spectralIndices(lsat, blue = 1, green=2, redEdge1=1, redEdge2=2, redEdge3=3, red=3, nir=4, swir2=5, swir3=7, coefs = list(L=0.4,s=0.3,swir2ccc=30,swir2coc=140), scaleFactor=255), "SpatRaster") - expect_equal(nlyr(sp), length(.IDXdb)) })