Smooths a 2D field using Discrete Cosine Transform or SVD

smooth2d.R

#' Smooths a 2D field
#' 
#' @param x,y Vector of x and y coordinates
#' @param value Vector of values
#' @param kx,ky Proportion of components to keep in the x and 
#' y direction respectively. Lower values increased the smoothness.
#' 
#' @examples 
#' library(ggplot2)
#' # Creates a noisy version of the volcano dataset and applies the smooth
#' volcano <- datasets::volcano |> 
#'   reshape2::melt(value.name = "original") |> 
#'   transform(noisy = original + 1.5*rnorm(length(original))) |> 
#'   transform(smooth_dct = smooth_dct(Var2, Var1, noisy, kx = 0.4)) |> 
#'   transform(smooth_svd = smooth_svd(Var2, Var1, noisy, variance_lost = 0.005)) |> 
#'   reshape2::melt(id.vars = c("Var1", "Var2"))
#' 
#' ggplot(volcano, aes(Var1, Var2)) +
#'   geom_contour(aes(z = value, color = ..level..)) +
#'   scale_color_viridis_c() +
#'   coord_equal() +
#'   facet_wrap(~variable, ncol = 2)
smooth_dct <- function (x, y, value, kx = 1, ky = kx) {
  force(ky)  # Evaluate ky now, because kx will change later.
  m <- mean(value)
  value <- value - m
  
  o <- order(y, x)
  # Transform into a matrix
  matrix <- matrix(value[o], nrow = length(unique(x)),  byrow = FALSE)
  
  # Compute the Discreete Cosine Transform
  # We could also use fft, but it suffers from edge effects
  dct <- gsignal::dct2(matrix)

  # Define the components to keep 
  # Remmeber that the FFT is symmetrical
  kx <- c(0, seq_len(floor(nrow(dct)/2 * kx)))
  kx <- c(kx + 1, nrow(dct) - kx[kx != 0] + 1)
  
  ky <- c(0, seq_len(floor(ncol(dct)/2 * ky)))
  ky <- c(ky + 1, ncol(dct) - ky[ky != 0] + 1)
  
  # Replace with zeros and invert
  dct[, -ky] <- 0
  dct[-kx, ] <- 0
  c(gsignal::idct2(dct))[order(o)] + m
}

#' @param  variance_lost Maximum percentage of variance lost after smoothing.
smooth_svd <- function(x, y, value, variance_lost = 0.01) {
  m <- mean(value)
  value <- value - m
  o <- order(y, x)
  # Transform into a matrix
  matrix <- matrix(value[o], nrow = length(unique(x)),  byrow = FALSE)
  total_variance <- norm(abs(matrix), type = "F")
  
  svd <- svd(matrix)
  
  variance_accumulated <- cumsum(svd$d^2/total_variance^2)
  variance_kept <- 1 - variance_lost
  
  keep <- seq_len(which(variance_accumulated - variance_kept > 0)[1])
  
  smooth <- svd$u[, keep] %*% diag(svd$d[keep], nrow = length(keep)) %*% t(svd$v[, keep])
  smooth <- smooth + m
  c(smooth)[order(o)]
}