Styjun

I have a matrix M (say r rows and c columns) and I would like to get the "weighted" sum for each matrix element based on it's neighbors and create a new matrix M2. The word neighbor could be within a radius of 1 (which is often called the Moore neighborhood in Cellular Automata theory) or the radius could be a different than 1, say, 2, 3, etc.

For a particular cell in the matrix M, say somewhere in the middle. Let's say position (i,j); then the (i,j)th cell has "eight" neighbors given by,

(i-1, j-1), (i-1, j), (i-1, j+1), (i, j-1), (i, j+1), (i+1, j-1), (i+1, j), (i+1, j+1).

I want to create a matrix M2 that calculates the "weighted" sum of the (i,j)th cell plus its eight neighbors. The weighting is done based on the Euclidean distance between cells. So for example,

exp(-sqrt(2))*M[i-1,j-1] + exp(-1)*M[i-1,j] + exp(-sqrt(2))*M[i-1,j+1] + exp(-1)*M[i,j-1] + M[i,j] + exp(-1)*M[i,j+1] + exp(-sqrt(2))*M[i+1,j-1] + exp(-1)*M[i+1,j] + exp(-sqrt(2))*M[i+1,j+1]

The same idea is repeated for all cells (cells along the boundaries need to be treated specially since they don't necessarily have eight neighboring cells). The above idea is for radius 1 but the code I am trying to develop needs to be generic for any radius.

r <- 4
c <- 4

n <- r*c

(M <- matrix(1:n, r, c))

addresses <- expand.grid(x = 1:r, y = 1:c)

#Got this code in the same forum

z <- rbind(c(-1,0,1,-1,1,-1,0,1),c(-1,-1,-1,0,0,1,1,1))

get.neighbors <- function(rw) 
 # Convert to absolute addresses 
 z2 <- t(z + unlist(rw))
 # Choose those with indices within mat 
 b.good <- rowSums(z2 > 0)==2 & z2[,1] <= nrow(M) & z2[,2] <= ncol(M)
 M[z2[b.good,]]


apply(addresses,1 , get.neighbors) # Returns a list with neighbors

M

Essentially, M2 for radius = 1 must be the "weighted" sum of the current cell plus the neighbors. The current current cell always gets a weight of 1.

M = [ 1 5 9 13
 2 6 10 14
 3 7 11 15
 4 8 12 16]

M2 = [ 5.033 13.803 .... ....
 .... .... .... ....
 .... .... .... ....
 .... .... .... ....]

How do I go about getting matrix M2 in R? What about if radius for more than 1? I would like the weighting to happening inside two for loops so I can use the calculated weighted sum of the [i,j] cell further in the code closing the two for loops.

I think the following does the weighted sum that you want. I'll be finding the neighbours in a similar way @r2evans did.

wtd_nbrs_sum <- function(input_matrix,
 radius,
 weight_matrix)

 temp_1 <- matrix(data = 0,
 nrow = nrow(x = input_matrix),
 ncol = radius)
 temp_2 <- matrix(data = 0,
 nrow = radius,
 ncol = ((2 * radius) + ncol(x = input_matrix)))
 input_matrix_modified <- rbind(temp_2,
 cbind(temp_1, input_matrix, temp_1),
 temp_2)
 output_matrix <- matrix(nrow = nrow(x = input_matrix),
 ncol = ncol(x = input_matrix))
 for(i in seq_len(length.out = nrow(x = input_matrix)))
 
 for(j in seq_len(length.out = nrow(x = input_matrix)))
 
 row_min <- (radius + (i - radius))
 row_max <- (radius + (i + radius))
 column_min <- (radius + (j - radius))
 column_max <- (radius + (j + radius))
 neighbours <- input_matrix_modified[(row_min:row_max), (column_min:column_max)]
 weighted_sum <- sum(neighbours * weight_matrix)
 output_matrix[i, j] <- weighted_sum
 
 
 return(output_matrix)


r <- 4
c <- 4
n <- r*c
M <- matrix(data = 1:n,
 nrow = r,
 ncol = c)
R <- 1
wts <- matrix(data = c(exp(x = -sqrt(x = 2)), exp(x = -1), exp(x = -sqrt(x = 2)), exp(x = -1), 1, exp(x = -1), exp(x = -sqrt(x = 2)), exp(x = -1), exp(x = -sqrt(x = 2))),
 nrow = 3,
 ncol = 3)

wtd_nbrs_sum(input_matrix = M,
 radius = R,
 weight_matrix = wts)
#> [,1] [,2] [,3] [,4]
#> [1,] 5.033856 13.80347 24.16296 23.89239
#> [2,] 8.596195 20.66391 34.43985 32.84175
#> [3,] 11.186067 24.10789 37.88383 35.43163
#> [4,] 9.748491 19.86486 30.22435 28.60703

^{Created on 2019-03-24 by the reprex package (v0.2.1)}

Hope this helps.

Edited to include weighted-sums.

There might be a really neat trick to doing it, but the most straight-forward (and maintainable) way is probably a simple two-for-loop implementation.

M1 <- matrix(1:16, nr=4)
M1
# [,1] [,2] [,3] [,4]
# [1,] 1 5 9 13
# [2,] 2 6 10 14
# [3,] 3 7 11 15
# [4,] 4 8 12 16

The code:

get_neighbors <- function(M, radius = 1) 
 M2 <- M
 M2[] <- 0
 nr <- nrow(M)
 nc <- ncol(M)
 eg <- expand.grid((-radius):radius, (-radius):radius)
 eg$wt <- exp(-sqrt(abs(eg[,1]) + abs(eg[,2])))
 for (R in seq_len(nr)) 
 for (C in seq_len(nc)) 
 ind <- cbind(R + eg[,1], C + eg[,2], eg[,3])
 ind <- ind[ 0 < ind[,1] & ind[,1] <= nr &
 0 < ind[,2] & ind[,2] <= nc,, drop = FALSE ]
 M2[R,C] <- sum(M[ind[,1:2, drop=FALSE]] * ind[,3])
 
 
 M2


get_neighbors(M1, 1)
# [,1] [,2] [,3] [,4]
# [1,] 5.033856 13.80347 24.16296 23.89239
# [2,] 8.596195 20.66391 34.43985 32.84175
# [3,] 11.186067 24.10789 37.88383 35.43163
# [4,] 9.748491 19.86486 30.22435 28.60703

The same thing, with radius 2:

get_neighbors(M1, 2)
# [,1] [,2] [,3] [,4]
# [1,] 12.44761 25.64963 31.73247 32.70974
# [2,] 18.57765 35.96237 43.33862 43.51911
# [3,] 20.09836 37.80643 45.18268 45.03982
# [4,] 17.51314 31.88500 37.96784 37.77527

And a simple test, if radius 0 is used, then M1 and M2 should be identical (they are).

Note: this generally performs just fine in base R, with no fancy use of apply or its cousins. Since this is a really straight-forward heuristic, it could easily be implemented with Rcpp to be significantly faster.