Applying a list of functions to a multiband raster using `terra::app()`

I have been processing some Sentinel 2 data using the "terra" and "RStoolbox" packages of R. I was trying to color balance two tiles of Sentinel 2 imagery using the histMatch() function, but got some very unexpected results...the maximum values for the bands histogram matched image were sometime less than the minimum values for the corresponding band in the reference image. I reported this behavior as a bug on GitHub as issue #117 on the RStoolbox repository. I dug into the source code and was able to isolate the issue. Essentially, histMatch() finds a function for each band in the reference raster that finds the raw value corresponding to a particular quantile in a cumulative distribution function. It then creates a list of functions and attempts to apply them sequentially to each corresponding band in the raster to be process. I pasted in the snippet of code from the raw histMatch() function below.

  totalFun <- function(xvals, f = layerFun) {
        if (is.vector(xvals)) 
            xvals <- as.matrix(xvals)
        app <- lapply(1:ncol(xvals), function(i) {
            f[[i]](xvals[, i])
        })
        do.call("cbind", app)
    }
    if (returnFunctions) {
        names(layerFun) <- names(x)
        return(layerFun)
    }
    .vMessage("Apply histogram match functions")
    out <- terra::app(x, fun = totalFun, ...)

Basically the problem is that terra::app() won't follow the order the functions are passed in. I am attempt to create a fix to post a proposed solution to my GitHub bug report. I have made a simpler problem that isolates the issue of passing multiple functions to terra::app() as a reproducible to demonstrate the desired output and the desire approach.

####Install and/or load terra####
if(!require(terra)){
  install.packages(terra)
}


####Example Data Setup####
set.seed(23)#For reproducibility

##Creating a Fake multiband raster##
r1a<-rast(xmin=-123.1, xmax=-122.1, ymin=45, ymax=45.5, nrows=100, ncols=100, crs="epsg:4326")
r1b<-rast(xmin=-123.1, xmax=-122.1, ymin=45, ymax=45.5, nrows=100, ncols=100, crs="epsg:4326")
r1c<-rast(xmin=-123.1, xmax=-122.1, ymin=45, ymax=45.5, nrows=100, ncols=100, crs="epsg:4326")

values(r1a)<-rnorm(10000, 200, 35)
values(r1b)<-runif(10000, 500, 10000)
values(r1c)<-rnorm(10000, 600, 150)

r1<-c(r1a, r1b, r1c)

##Creating a list of functions
fxn1<-function(x){log(x+1)}
fxn2<-function(x){(x)^0.5}
fxn3<-function(x){x-25}

funtime<-list(fxn1, fxn2, fxn3)

####Applying each function to its corresponding band in the fake raster####
##Desired behavior##
out<-list()
for(i in 1:nlyr(r1)){
  out[[i]]<-app(r1[[i]], funtime[[i]])
}

do.call(c, out)

####Desired approach####
##Will throw error##
totalfun<-function(x, f=funtime){
  lapply(1:length(x), function(i){
  out<-f[[i]](x[[i]])
})
  do.call(c, out)
}

result<-app(r1, fun=totalfun)

If anyone can think of a workaround to this problem I would greatly appreciate it, and it may help others attempting to histogram match other multiband imagery. Thanks!

Your example output:

x <- do.call(c, out)

Can be obtained with

f1 <- function(x) {
    cbind(fxn1(x[,1]), fxn2(x[,2]), fxn3(x[,3]))
}
y <- app(r1, f1)

all(values(x) == values(y))
#[1] TRUE

And with

ff <- list(fxn1, fxn2, fxn3)
f2 <- function(x) {
    sapply(1:ncol(x), \(i) ff[[i]](x[,i]))
}
z <- app(r1, f2)

all(values(x) == values(z))
#[1] TRUE