R: Apply FUN to kxk subsections of array

Question

The language is R.

I have an nxm matrix, and I'd like to partition it into 3x3 sections and calculate the mean (or any function) within each. (If there's a leftover bit that isn't 3x3 then use just what's left).

I'm sure there's an apply-ish way to do this -- it's on the tip of my tongue -- but my brain is currently failing me. I suppose it's a bit like a moving window question except I want non-overlapping windows (so it's easier).

Can anyone think of an inbuilt function that does this? Or a vectorised way?

Here's my loopy version:

winSize <- 3
mat <- matrix(runif(6*11),nrow=6,ncol=11)
nr <- nrow(mat)
nc <- ncol(mat)
outMat <- matrix(NA,nrow=ceiling(nr/winSize),
                    ncol=ceiling(nc/winSize))
FUN <- mean
for ( i in seq(1,nr,by=winSize) ) {
    for ( j in seq(1,nc,by=winSize) ) {
        # work out mean in 3x3 window, fancy footwork
        #  with pmin just to make sure we don't go out of bounds
        outMat[ ceiling(i/winSize), ceiling(j/winSize) ] <-
               FUN(mat[ pmin(i-1 + 1:winSize,nr), pmin(j-1 + 1:winSize,nc)])
    }
}

cheers.

Answer 1

You can use row and col to extract the row and column numbers, and then compute the coordinates of each block.

tapply( 
  mat, 
  list( floor((row(mat)-1)/winSize), floor((col(mat)-1)/winSize) ), 
  mean 
)

Edit: This can be generalized to higher-dimensional arrays, by replacing row and col with the following function.

a <- function( m, k ) {
  stopifnot( "array" %in% class(m) || "matrix" %in% class(m) )
  stopifnot( k == floor(k) )
  stopifnot( k > 0 )
  n <- length(dim(m))
  stopifnot( k <= n )
  i <- rep(
    1:dim(m)[k],
    each  = prod(dim(m)[ 1:n < k ]),
    times = prod(dim(m)[ 1:n > k ])
  )  
  array(i, dim=dim(m))
}

# A few tests
m <- array(NA, dim=c(2,3))
all( row(m) == a(m,1) )
all( col(m) == a(m,2) )
# In dimension 3, it can be done manually:
m <- array(NA, dim=c(2,3,5))
all( a(m,1) == array( rep(1:dim(m)[1], times=prod(dim(m)[2:3])), dim=dim(m) ) )
all( a(m,2) == array( rep(1:dim(m)[2], each=dim(m)[1], times=dim(m)[3]), dim=dim(m) ) )
all( a(m,3) == array( rep(1:dim(m)[3], each=prod(dim(m)[-3])), dim=dim(m) ) )