How to generate a discrete 2D Gaussian smoothing kernel using R

Question

I'm wondering whether R has a function for this, or how to code, a 2D discrete version of a Gaussian smoothing kernel. The goal is to use this to smooth a 2D matrix of values.

Answer 1

I wrote a function to solve my problem. Here is the function that I've called kernelsmooth. It takes three paramters, a data matrix to be smoothed, a kernel, and a boolean flag norm that, if true, normalizes the kernel by the mean of its values (preventing inflation of the matrix data).

kernelsmooth <- function(x, kern, norm=TRUE) {
  # how many rows/cols of zeroes are used to pad.
  width <- dim(kern)[1]
  pad <- floor(width / 2)

  # record the width and height the input data matrix
  x_w <- ncol(x)
  x_h <- nrow(x)

  # Are we normalizing the kernel?
  if (norm == TRUE) {
    k <- kern / sum(abs(kern))
  } else {
    k <- kern
  }

  # pad all around the matrix an equal width of zeros
  x_pad <- t(padzeros(data=x, nzeros=pad, side="both"))
  x_pad <- t(padzeros(data=x_pad, nzeros=pad, side="both"))

  # Pre-allocate the final (smoothed) data matrix
  s <- matrix(0, nrow = x_h, ncol = x_w)

  # Pre-allocate a temporary matrix for the iterative calculations
  temp <- matrix(0, width, width)

  # Loop through the data to apply the kernel.
  for (col in 1:x_w ) {
    for (row in 1:x_h ) {
      temp <- x_pad[row:(row + width - 1), col:(col + width - 1)]
      s[row,col] <-  sum(k * temp)
    }
  }

  # return the smoothed data
  return(s)
}

Some examples:

1. Do not smooth a matrix (apply an identity kernel analogous to an identity matrix).

   kernelsmooth(x = matrix(2,5,5), kern = matrix(c(0,0,0, 0,1,0, 0,0,0) , 3, 3))

Output:

     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    2    2    2    2    2    2
[2,]    2    2    2    2    2    2
[3,]    2    2    2    2    2    2
[4,]    2    2    2    2    2    2
[5,]    2    2    2    2    2    2
[6,]    2    2    2    2    2    2

1. Average smoothing (each element is smoothed by the average if itself and its neighbours):

   kernelsmooth(x = matrix(4,5,5), kern = matrix(1,3,3))

Output:

         [,1]     [,2]     [,3]     [,4]     [,5]
[1,] 1.777778 2.666667 2.666667 2.666667 1.777778
[2,] 2.666667 4.000000 4.000000 4.000000 2.666667
[3,] 2.666667 4.000000 4.000000 4.000000 2.666667
[4,] 2.666667 4.000000 4.000000 4.000000 2.666667
[5,] 1.777778 2.666667 2.666667 2.666667 1.777778

1. Gaussian smoothing (isotropic smoothing):

   gauss <- ( c(1,4,7,4,1) %*% t(c(1,4,7,4,1)) ) kernelsmooth(x = matrix(5,8,8), kern = gauss)

Output:

         [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]
[1,] 2.491349 3.321799 3.529412 3.529412 3.529412 3.529412 3.321799 2.491349
[2,] 3.321799 4.429066 4.705882 4.705882 4.705882 4.705882 4.429066 3.321799
[3,] 3.529412 4.705882 5.000000 5.000000 5.000000 5.000000 4.705882 3.529412
[4,] 3.529412 4.705882 5.000000 5.000000 5.000000 5.000000 4.705882 3.529412
[5,] 3.529412 4.705882 5.000000 5.000000 5.000000 5.000000 4.705882 3.529412
[6,] 3.529412 4.705882 5.000000 5.000000 5.000000 5.000000 4.705882 3.529412
[7,] 3.321799 4.429066 4.705882 4.705882 4.705882 4.705882 4.429066 3.321799
[8,] 2.491349 3.321799 3.529412 3.529412 3.529412 3.529412 3.321799 2.491349

Answer 2

The sharp drop-off at the edges of your matrices, e.g. from 3.529412 to 0, is bad news. You need to enlarge the matrix to such a size that the edge values are as close to zero as possible. You can check the validity of this statement by computing the 2D Fourier transform of each of your matrices, and comparing this with the same transform done on the larger versions. You will also see much better results when applied to your data.