R/coregister.R

In mmontesinosanmartin/rsatExtra: What the package does (short line)

# provisional: Kudos to RSToolbox
# TODO: implementation in c++
# TODO: add rotation in addition to x-y shifts
# TODO: work with rtois
coregister <- function(s1, s2, xy.shift = 3, shiftInc = 1, nSamples = 1e5, reportStats = TRUE, verbose, nBins = 100, ...) {
  
  if(!missing("verbose")) .initVerbose(verbose)
  if(!compareCRS(s2,  s1)) stop("Projection must be the same for s2 and s1")
  nSamples <- min(nSamples, ncell(s1))
  
  if(inherits(xy.shift, "matrix") && ncol(xy.shift)  == 2) { 
    xy.shifts <- xy.shift # * res(s1) 
  } else {
    xy.shift <-  seq(0, xy.shift, shiftInc)
    xy.shift <- c(-rev(xy.shift), xy.shift[-1])
    # xy.shifts <- expand.grid(xy.shift * res(s1)[1], xy.shift * res(s1)[2])
    xy.shifts <- expand.grid(xy.shift , xy.shift)
  } 
  names(xy.shifts) <- c("x", "y")        
  
  ran   <- apply(xy.shifts, 2, range)
  minex <- extent(shift(s1, ran[1,1], ran[1,2]))
  maxex <- extent(shift(s1, ran[2,1], ran[2,2]))   
  
  XYs1s <- sampleRandom(s2, size = nSamples, ext = .getExtentOverlap(minex, maxex)*0.9, xy = TRUE)
  xy <- XYs1s[,c(1,2)]
  me <- XYs1s[,-c(1,2)]     
  mmin <- min(minValue(s2))
  mmax <- max(maxValue(s2))
  smin <- min(minValue(s1))
  smax <- max(maxValue(s1))
  
  mbreax <- seq(mmin, mmax, by = (mmax - mmin)/nBins)
  sbreax <- seq(smin, smax, by = (smax - smin)/nBins)
  me       <- cut(me, breaks = mbreax, labels = FALSE, include.lowest = TRUE)
  
  nsl <- nlayers(s1)
  nml <- nlayers(s2)
  if(nsl !=  nml)  stop("Currently s1 and s2 must have the same number of layers")         
  
  shiftPts <- function(o, x, y) {
    o[,"x"] <- o[,"x"] + x
    o[,"y"] <- o[,"y"] + y
    o
  }
  
  spts <- .parXapply(X = 1:nrow(xy.shifts), XFUN = "lapply", FUN = function(i){
    xt <- shiftPts(xy, x = -xy.shifts[i,1], y = -xy.shifts[i,2])
    cellFromXY(s1, xt)
  }, envir=environment())
  
  ucells <- sort(unique(unlist(spts)))
  lut <- as.matrix(s1[ucells])
  rownames(lut) <- ucells
  spts <- lapply(spts, as.character)
  
  sh <- .parXapply(X = 1:nrow(xy.shifts), XFUN = "lapply", FUN = function(i){
    se <- lut[spts[[i]], ]
    se  <- cut(se, breaks = sbreax, labels = FALSE, include.lowest = TRUE)  
    
    pab  <- table(me,se)
    pab  <- pab/sum(pab)
    
    pa   <- colSums(pab)
    pb   <- rowSums(pab)
    
    pabx <- pab[pab>0]
    pb   <- pb[pb>0] 
    pa   <- pa[pa>0]
    
    hab  <- sum(-pabx * log(pabx))
    ha   <- sum(-pa * log(pa))
    hb   <- sum(-pb * log(pb))
    mi   <- ha + hb -hab
    list(mi = mi, joint = pab)
    
  }, envir = environment() )
  
  mi <- vapply(sh,"[[", i = 1, numeric(1)) 
  if(reportStats){
    jh <- lapply(sh, function(x) matrix(as.vector(x$joint), nrow=nrow(x$joint), ncol=ncol(x$joint)))   
    names(jh) <- paste(xy.shifts[,"x"], xy.shifts[,"y"], sep = "/")
  }
  moveIt <- xy.shifts[which.max(mi),]
  moved <- shift(s1, dx = moveIt[1], dy = moveIt[2])
  return(list(MI = data.frame(xy.shifts, mi=mi), jointHist = jh, optimal = moveIt, img = moved))
}

.getExtentOverlap <- function (...) 
{
  el <- list(...)
  if (any(!vapply(el, inherits, what = "Extent", logical(1)))) 
    stop("You can only supply Extent objects to getExtentOverlap")
  em <- do.call("rbind", lapply(el, as.vector))
  extent(c(max(em[, 1]), min(em[, 2]), max(em[, 3]), min(em[, 
                                                            4])))
}

.parXapply <- function (X, XFUN, MARGIN, FUN, envir, ...) 
{
  call <- quote(f(cl = cl, X = X, FUN = FUN, MARGIN = MARGIN, 
                  ...))
  if (isTRUE(getOption("rasterCluster"))) {
    cl <- getCluster()
    on.exit(returnCluster())
    f <- c(lapply = parLapply, sapply = parSapply, apply = parApply)[[XFUN]]
    if (!is.primitive(FUN)) {
      g <- findGlobals(FUN)
      gg <- lapply(g, get, envir = envir)
      names(gg) <- g
    }
    else {
      gg <- NULL
    }
    l <- c(list(...), gg)
    clusterExport(cl = cl, names(l), envir = envir)
    if (XFUN == "lapply") 
      names(call)[names(call) == "FUN"] <- "fun"
  }
  else {
    f <- get(XFUN)
    call[["cl"]] <- NULL
  }
  if (XFUN != "apply") 
    call[["MARGIN"]] <- NULL
  eval(call)
}