R: Применить FUN к kxk подразделам массива

Вы можете использовать row и col для извлечения номеров строк и столбцов, а затем вычислить координаты каждого блока.

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

Изменить: это можно обобщить на массивы более высокой размерности, заменив row и col следующей функцией.

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) ) )

Просто хочу обобщить различные методы для этого.

Во-первых, решение @VincentZoonekynd. Это очень общий способ — он позволяет мне применять любую функцию к моей матрице. Однако это немного медленно, потому что я применяю их к матрицам порядка ~ 5000x1000x3 и хочу вернуть изображение (5000/kernelSize) x (1000/kernelSize) x 3.

Во-первых, создайте матрицу для тестирования (я сделал ее меньше, чтобы не убить свой компьютер при тестировании различных методов):

sz <- c(1000,300,3)
img <- array(runif(prod(sz)),dim=sz)
kernelSize <- 3
outSz <- c(ceiling(sz[1:2]/kernelSize),3)
FUN <- mean

Метод 0: цикл

############
# METHOD 0 #
############
# Loopy. base standard.
t0 <- system.time({
out0 <- array(NA,dim=outSz)
for ( i in seq(1,sz[1],by=kernelSize) ) {
    for ( j in seq(1,sz[2],by=kernelSize) ) {
        for ( c in 1:sz[3] ) {
        # work out mean in 3x3 window, fancy footwork
        #  with pmin just to make sure we don't go out of bounds
        out0[ ceiling(i/kernelSize), ceiling(j/kernelSize),c ] <-
               FUN(img[ pmin(i-1 + 1:kernelSize,sz[1]), 
                        pmin(j-1 + 1:kernelSize,sz[2]),
                        c]) 
        }
    }
}})

Способ 1: нажмите

############
# METHOD 1 #
############
# @Vincent Zoonekynd.
# I can apply *any* function I want. how awesome!
# NOTE: I just realised that there is a slice.index(img,i)
#       is the same as his a(img,i) function.
t1 <- system.time({
out1 <- tapply(
         img,
         list( floor((slice.index(img,1)-1)/kernelSize), 
               floor((slice.index(img,2)-1)/kernelSize),
               slice.index(img,3) ),
         FUN )
})

cat('METHOD 0:',t0['elapsed'],'\n')
cat('METHOD 1:',t1['elapsed'],'\n')
cat(all(out0==out1),'\n')

Это дает:

METHOD 0: 13.549 
METHOD 1: 19.415 
TRUE

Что немного медленно, учитывая, что я хотел бы применить это к большим матрицам img.

Что меня удивило (сначала), так это то, что МЕТОД 0 (циклы) был быстрее, чем МЕТОД 1 (tapply).

Однако затем я вспомнил, что tapply имеет репутацию не намного быстрее, чем явный цикл (почему? Я помню, что где-то читал об этом... код функции выглядит так, как будто он все равно может выполнять цикл for, в отличие от вызова внешнего код).

У меня также есть общее ощущение, что vapply и sapply являются быстрыми версиями apply (опять же, не уверен, что это окончательно верно, но я определенно нашел это).

Способ 2: ваппли

Итак, я попытался переписать свою зацикленную версию, используя vapply. (Возможно, есть лучший способ справиться с 3-м измерением, но да ладно...). Это в основном генерирует большой список координат в img. Координаты дают (i,j) угла каждого kernelSize*kernelSize квадрата.

Затем vapply перебирает их все и вычисляет среднее значение.

##########
# METHOD 2 
##########
# use 'vapply'
t2 <- system.time({
is <- seq(1,sz[1],by=kernelSize)
js <- seq(1,sz[2],by=kernelSize)
# generate a (nrow*nsize) x 2 array with
# all (i,j) combinations for corners of
# kernelSize*kernelSize squares.
# Do it column-major so we can reshape after.
coords <- cbind( rep.int(is,length(js)), rep(js,each=length(is)) ) 
out2 <- array(NA,dim=outSz)
for ( c in 1:sz[3] ) { 
    out2[,,c] <- array(
    vapply( 1:nrow(coords), function(i) {
          FUN(img[coords[i,1]:pmin(sz[1],coords[i,1]+kernelSize-1),
                   coords[i,2]:pmin(sz[2],coords[i,2]+kernelSize-1),
                   c])
            }, 0 ),
                dim=outSz[1:2] ) 
}})
cat('METHOD 2:',t2['elapsed'],'\n')
cat(all(out0==out2),'\n')

Это дает:

METHOD 2: 12.627 
TRUE

Таким образом, это немного быстрее, чем цикл, для использования vapply (я чувствую, что не получаю столько от vapply, как мог бы, хотя... как будто я не использую его правильно).

Способ 3: фильтр

Это по-прежнему недостаточно быстро, поэтому я включил информацию о том, что мне требуется только среднее значение в каждом окне, и это в основном свертка [ 1/3 1/3 1/3 ] с матрицей в каждом измерении.

Это теряет общую применимость применения произвольного FUN, но взамен дает большое ускорение.

По сути, я создаю ядро ​​[1/3, 1/3, 1/3] и дважды сворачиваю его с помощью img, один раз в направлении x и один раз в направлении y. Затем я извлекаю только каждое третье значение (так как мне нужны неперекрывающиеся окна).

Мне это кажется немного расточительным, поскольку я вычисляю среднее значение для каждого окна 3x3 в моей исходной матрице, а не только для неперекрывающихся окон, но я не знаю как сказать R не вычислять те значения, которые я все равно собираюсь выбросить.

Однако вы должны проявлять осторожность на границах - скажем, остается только участок 2x2, тогда среднее значение превышает 4 вместо 9 значений. Мой текущий код не обрабатывает это, но я не возражаю, если это просто граница, потому что я делаю субдискретизацию только для целей отображения.

(Хотя было бы неплохо исправить эту последнюю вещь...)

##########
# METHOD 3 
##########
# Convolve using `filter`,
# since the mean in a window is just a 
# convolution.
t3 <- system.time({
is <- pmin(seq(1,sz[1],by=kernelSize) + floor(kernelSize/2),sz[1]-1)
js <- pmin(seq(1,sz[2],by=kernelSize) + floor(kernelSize/2),sz[2]-1)
out3 <- array(NA,dim=outSz)
for ( c in 1:3 ) {
    out3[,,c] <- (t(filter(
                    t(filter(img[,,c],rep(1,kernelSize))),
                    rep(1,kernelSize))))[is,js]
}
out3 <- out3/(kernelSize*kernelSize)
})
cat('METHOD 3:',t3['elapsed'],'\n')
cat(sum(out0!=out3),'\n')

Это возвращает:

METHOD 3: 1.593 
300

Таким образом, этот метод намного самый быстрый, и ошибка только в последнем столбце out3 (один раз на канал), поскольку (я думаю) существуют граничные условия.