I am trying to perform 2D matrix processing (convolution) by applying a kernel/filter to a large matrix. There is a built-in function convolve that can perform convolution. It offers three different types of convolution, namely circular, open and filter. It looks like filter is what I am looking for, but it isn't.
Let's consider two matrices:
# "large" matrix
(largemx <- matrix(1:16, nrow=4, byrow = T))
[,1] [,2] [,3] [,4]
[1,] 1 2 3 4
[2,] 5 6 7 8
[3,] 9 10 11 12
[4,] 13 14 15 16
# kernel
(kernel <- rbind(1,c(1, 0, 1), 1))
[,1] [,2] [,3]
[1,] 1 1 1
[2,] 1 0 1
[3,] 1 1 1
This is what I am getting with the filter argument (a vector of 8 elements):
convolve(largemx, kernel, type="f")
[1] 61 63 65 67 69 71 73 75
With type="open" the output is a vector of 17 elements:
convolve(largemx, kernel, type="o")
[1] 1 6 15 28 29 31 37 47 61 63 65 67 69 71 73 75 72 65 54 39 40 36 28 16
The circular type convolution requires the kernel to be the same size as the matrix and the output is a matrix too:
(cbind(rbind(kernel,0),0) -> kernelc)
[,1] [,2] [,3] [,4]
[1,] 1 1 1 0
[2,] 1 0 1 0
[3,] 1 1 1 0
[4,] 0 0 0 0
convolve(largemx, kernelc, type="c")
[,1] [,2] [,3] [,4]
[1,] 48 56 52 52
[2,] 80 88 84 84
[3,] 64 72 68 68
[4,] 64 72 68 68
Below follows the desired output. I would very much appreciate if someone could clarify whether any of those types of convolution could be used for the task and maybe explain how the calculations are performed.
[,1] [,2] [,3] [,4]
[1,] 13 22 27 18
[2,] 28 48 56 37
[3,] 48 80 88 57
[4,] 33 58 63 38
Inspired by this question on CGCC
I have found an easy and compact solution how to apply convolve for this case. This is a bit hacky but works for my particular task.
This is the input matrix once again:
largemx <- matrix(1:16, nrow=4, byrow = T)
apply the filter c(1, 1, 1) once to the rows padded with zeroes:
(convoluted_rows <- apply(largemx, 1, \(x)convolve(c(0, x, 0), c(1, 1, 1), type="f")))
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0 3 11 19 27 0
[2,] 0 6 18 30 42 0
[3,] 0 9 21 33 45 0
[4,] 0 7 15 23 31 0
And then again to the result:
(convoluted_twice <- apply(convoluted_rows, 1, \(x)convolve((0, x, 0), c(1, 1, 1), type="f")))
[,1] [,2] [,3] [,4]
[1,] 14 24 30 22
[2,] 33 54 63 45
[3,] 57 90 99 69
[4,] 46 72 78 54
It only remains to subtract the original matrix:
(convoluted_twice-largemx)
[,1] [,2] [,3] [,4]
[1,] 13 22 27 18
[2,] 28 48 56 37
[3,] 48 80 88 57
[4,] 33 58 63 38
How this works? (in case you were wondering)
As it is written in the description to convolve fuction, "Use the Fast Fourier Transform to compute the several kinds of convolutions of two sequences." So it is designed to operate on 1D and not on 2D objects. When you are feeding 2D arrays to convolve, they get flattened, which is unwanted.
In my solution I have started by convolving the rows (by the way, it will work with the columns the same way, but you need to transpose the result) with the kernel c(1, 1, 1). It sums up the rows by the groups of consequent three numbers: Σ[0, 1, 2] = 3, Σ[1, 2, 3] = 6, Σ[2, 3, 4] = 9 etc. Note that apply gives a transposed output, but since it is used twice, that doesn't matter.
In the second run the same kernel is applied to the result. Now the columns are summed up by groups of three. In the end this produces a matrix where each element is a sum of the corresponding 3x3 sumbatrix of pad (see the jay.sf's answer). Because the host number of its cell in the matrix should not be included, we simply subtract the original largemx from the result.