#
# blur.R
#
# apply Gaussian blur to an image
#
#    $Revision: 1.12 $   $Date: 2009/12/16 18:23:52 $
#
fillNA <- function(x, value=0) {
  stopifnot(is.im(x))
  v <- x$v
  v[is.na(v)] <- value
  x$v <- v
  return(x)
}
  
blur <- function(x, sigma=NULL, ..., normalise=FALSE, bleed=TRUE, varcov=NULL) {
  stopifnot(is.im(x))
  # determine smoothing kernel 
  sigma.given <- !is.null(sigma)
  varcov.given <- !is.null(varcov)
  if (sigma.given) {
    stopifnot(is.numeric(sigma))
    stopifnot(length(sigma) %in% c(1, 2))
    stopifnot(all(sigma > 0))
  }
  if (varcov.given)
    stopifnot(is.matrix(varcov) && nrow(varcov) == 2 && ncol(varcov) ==
              2)
  ngiven <- varcov.given + sigma.given
  switch(ngiven + 1,
         {
           sigma <- (1/8) * min(diff(x$xrange), diff(x$yrange))
         }, {
           if (sigma.given && length(sigma) == 2)
             varcov <- diag(sigma^2)
           if (!is.null(varcov))
             sigma <- NULL
         }, {
           stop(paste("Give only one of the arguments", sQuote("sigma"),
                      "and", sQuote("varcov")))
         })
  # replace NA's in image raster by zeroes 
  X <- fillNA(x, 0)
  # convolve with Gaussian
  Y <- second.moment.calc(X, sigma=sigma, varcov=varcov, what="smooth")
  # if no bleeding, we restrict data to the original boundary
  if(!bleed)
    Y$v[is.na(x$v)] <- NA
  # 
  if(!normalise)
    return(Y)
  # normalisation:
  # convert original image to window (0/1 image)
  Xone <- x
  isna <- is.na(x$v)
  Xone$v[isna] <- 0
  Xone$v[!isna] <- 1
  # convolve with Gaussian
  Ydenom <- second.moment.calc(Xone, sigma=sigma, ..., varcov=varcov, what="smooth")
  # normalise
  Z <- eval.im(Y/Ydenom)
  return(Z)
}
  
safelookup <- function(Z, X, factor=2, warn=TRUE) {
  # X is a ppp
  # evaluates Z[X], replacing any NA's by blur(Z)[X]
  Zvals <- Z[X, drop=FALSE]
  if(any(isna <- is.na(Zvals))) {
    # First pass - look up values at neighbouring pixels if valid
    XX <- X[isna]
    rc <- nearest.valid.pixel(XX$x, XX$y, Z)
    Zvals[isna] <- Z$v[cbind(rc$row, rc$col)]
  }
  if(any(isna <- is.na(Zvals))) {
    # Second pass - extrapolate
    XX <- X[isna]
    pixdiam <- sqrt(Z$xstep^2 + Z$ystep^2)
    # expand domain of Z 
    RX <- as.rectangle(X)
    RZ <- as.rectangle(Z)
    bb <- bounding.box(RX, RZ)
    big <- grow.rectangle(bb, 2 * pixdiam)
    Z <- rebound.im(Z, big)
    # now blur
    Zblur <- blur(Z, factor * pixdiam, bleed=TRUE, normalise=TRUE)
    Bvals <- Zblur[XX, drop=FALSE]
    if(any(is.na(Bvals))) 
      stop("Internal error: pixel values were NA, even after blurring")
    Zvals[isna] <- Bvals
    if(warn)
      warning(paste(sum(isna), "out of", X$n, "pixel values",
                    "were outside the pixel image domain",
                    "and were estimated by convolution"))
  }
  return(Zvals)
}