Last active
August 10, 2022 16:57
-
-
Save eliocamp/5850c629a032fbeae39931d236e82fc7 to your computer and use it in GitHub Desktop.
Smooths a 2D field using Discrete Cosine Transform or SVD
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#' 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)] | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment