R/pifMatch.R

In RStoolbox: Remote Sensing Data Analysis

#' Pseudo-Invariant Features based Image Matching
#' 
#' Match one scene to another based on linear regression of pseudo-invariant features (PIF).
#' 
#' @param img SpatRaster. Image to be adjusted.
#' @param ref SpatRaster. Reference image.
#' @param method Method to calculate pixel similarity. Options: euclidean distance ('ed'), spectral angle ('sam') or pearson correlation coefficient ('cor').
#' @param quantile Numeric. Threshold quantile used to identify PIFs
#' @param returnPifMap Logical. Return a binary raster map ot pixels which were identified as pesudo-invariant features.
#' @param returnSimMap Logical. Return the similarity map as well
#' @param returnModels Logical. Return the linear models along with the adjusted image.
#' @details 
#' The function consists of three main steps:
#' First, it calculates pixel-wise similarity between the two rasters and identifies pseudo-invariant pixels based on 
#' a similarity threshold. 
#' In the second step the values of the pseudo-invariant pixels are regressed against each other in a linear model for each layer.
#' Finally the linear models are applied to all pixels in the \code{img}, thereby matching it to the reference scene.
#' 
#' Pixel-wise similarity can be calculated using one of three methods: euclidean distance (\code{method = "ed"}), spectral angle (\code{"sam"}) or pearsons correlation coefficient (\code{"cor"}).
#' The threshold is defined as a similarity quantile. Setting \code{quantile=0.95} will select all pixels with a similarity above the 95\% quantile as pseudo-invariant features.
#' 
#' Model fitting is performed with simple linear models (\code{\link[stats]{lm}}); fitting one model per layer. 
#' @return 
#' Returns a List with the adjusted image and intermediate products (if requested). 
#' \itemize{
#'    \item \code{img}: the adjusted image
#'    \item \code{simMap}: pixel-wise similarity map (if \code{returnSimMap = TRUE})
#'    \item \code{pifMap}: binary map of pixels selected as pseudo-invariant features (if \code{returnPifMap = TRUE}) 
#'    \item \code{models}: list of linear models; one per layer (if \code{returnModels = TRUE})                          
#' }
#' @export 
#' @examples 
#' library(terra)
#'
#' 
#' ## Create fake example data
#' ## In practice this would be an image from another acquisition date
#' lsat_b <- log(lsat)
#' 
#' ## Run pifMatch and return similarity layer, invariant features mask and models
#' lsat_b_adj <- pifMatch(lsat_b, lsat, returnPifMap = TRUE,
#'                          returnSimMap = TRUE, returnModels = TRUE)
#' \donttest{
#' ## Pixelwise similarity
#' ggR(lsat_b_adj$simMap, geom_raster = TRUE)
#' 
#' ## Pesudo invariant feature mask 
#' ggR(lsat_b_adj$pifMap)
#' 
#' ## Histograms of changes
#' par(mfrow=c(1,3))
#' hist(lsat_b[[1]], main = "lsat_b")
#' hist(lsat[[1]], main = "reference")
#' hist(lsat_b_adj$img[[1]], main = "lsat_b adjusted")
#' 
#' ## Model summary for first band
#' summary(lsat_b_adj$models[[1]])
#' }
pifMatch <- function(img, ref, method = "cor", quantile = 0.95, returnPifMap = TRUE, returnSimMap = TRUE, returnModels = FALSE){
    img <- .toTerra(img)
	ref <- .toTerra(ref)

     if(nlyr(img)!=nlyr(ref) | nlyr(img) <= 1)
      stop("Both images need at least two corresponding bands and must have the same number of bands.", call.=FALSE)

    imgfull <- img
    ## Get joint extent
    if(!ext(img)==ext(ref)) {
        img <- crop(img, ref)
        ref <- crop(ref, img)
    }
    
    ## Calculate pixelwise similarity
    if(!method %in% c("ed", "sam", "cor"))
      stop("method must be one of 'ed', 'cor' or 'sam'", call. = FALSE)

    nmeth <- c(ed=1, sam=2, cor=3)[method]

    pifield_values <- pwSimilarityCpp(img[], ref[] ,nmeth)
    pifield <- img[[1]]
    values(pifield) <- pifield_values
    names(pifield) <- "layer"

    names(pifield) <- method

    ## Similarity quantile threshold
    thresh  <- quantile(values(pifield), p = quantile)
    pi      <- pifield > thresh

    names(pi) <- "pifMap"

    ## Fit linear models
    models <- lapply(1:nlyr(img), function(i){
        df <- data.frame(img = img[[i]][pi], ref = ref[[i]][pi])
        colnames(df) <- c("img", "ref")
        mod <- lm(ref ~ img, data = df)
    })

    ## Predict linear models for whole raster
    correct <- lapply(seq_along(models), function(i){
        corLayer <- imgfull[[i]]
        mod <- models[[i]]
        names(corLayer) <- "img"
        predict(corLayer, mod)
    })
    
    ## Assemble and return output
    names(models) <- names(img)
    names(correct) <- names(imgfull)
    
    out <- list( img = correct, simMap = pifield, pifMap = pi, models = models )
    ret <- c("img",
            if(returnSimMap) "simMap",
            if(returnPifMap) "pifMap",
            if(returnModels) "models"
    )
    out[!names(out) %in% ret] <- NULL
    return(out)
}