Блог

Я работаю с кубами изображений и хотел бы извлечь указанные отдельные пиксели из всех слоев (куб изображения похож на стек растров в том, что в нем есть N слоев, каждый из которых соответствует определенной длине волны; каждый слой представляет собой полное 2D-изображение). Я предложил два подхода:

 # sample 3-D image cube array
foo<- array(runif(512*512*16),dim=c(512,512,16))
# sample collection of x,y ordered pairs defining the pixels of interest
fooloc<- cbind(sample(1:512,256),sample(1:512,256))

pull1<-function(foo,fooloc) {
bar<-foo[fooloc[,1],fooloc[,2],]
sapply(1:(nrow(fooloc)),function(j)diag(bar[,,j]))
}

pull2 <- function(foo, fooloc) sapply(1:(nrow(fooloc)) ,function(j) foo[fooloc[j,1],fooloc[j,2],])
  

Второй примерно в 10 раз быстрее (на microbenchmark ). У меня просто такое ощущение, что я чего-то не хватает в возможностях [ оператора. Есть идеи?

РЕДАКТИРОВАТЬ: попытка raster идеи Джбаума:

 bfoo<-brick(foo)
pull3<-function(foo,fooloc) {
    #foo<-brick(foo)
    locfoo<- nrow(foo)*(fooloc[,1]-1)   fooloc[,2]
    extract(foo,locfoo) #default is all layers
    }
  

эталон для (фиксированной) pull2 функции:

 Unit: microseconds
                expr     min      lq   median       uq      max
 pull3(bfoo, fooloc) 924.523 932.221 1028.435 1041.691 1380.369
  pull2(foo, fooloc) 771.862 793.670  828.307  833.867  979.687
 neval
     5
     5
  

Если столбцы fooloc представляют местоположения строк (x) и столбцов (y), в которых вы хотели бы выполнить выборку (по всем слоям) foo , то способ, которым вы foo создали подмножество bar , для меня не имеет смысла. Выбранные строки foo , которые сохраняются в bar , — это строка и номер столбца первого интересующего местоположения. Выбранные столбцы, которые сохраняются bar , — это номер строки и столбца второго интересующего местоположения. Это происходит в обоих pull1 и pull2 .

Предполагая, что вы хотите выполнить детализацию по слоям в каждой точке (представленной координатой x, заданной в первом столбце, fooloc и координатой y, заданной во втором столбце fooloc ), тогда вы можете использовать матричное подмножество, посредством которого вы передаете матрицу со столбцами для каждого измерения (x, y, zв вашем случае), и строки, представляющие интересующие местоположения (т. Е., В основном, Ваш fooloc с дополнительным столбцом для номера среза):

 xyz <- cbind(fooloc[rep(1:256, each=16), ], rep(1:16, 256))

head(xyz, 17)

#       [,1] [,2] [,3]
#  [1,]  326  264    1
#  [2,]  326  264    2
#  [3,]  326  264    3
#  [4,]  326  264    4
#  [5,]  326  264    5
#  [6,]  326  264    6
#  [7,]  326  264    7
#  [8,]  326  264    8
#  [9,]  326  264    9
# [10,]  326  264   10
# [11,]  326  264   11
# [12,]  326  264   12
# [13,]  326  264   13
# [14,]  326  264   14
# [15,]  326  264   15
# [16,]  326  264   16
# [17,]  244  355    1


baz <- foo[xyz]

# Above, we use fooloc[rep(1:256, each=16), ] to repeat each location 
#  (each row of fooloc) 16 times, and bind a new column to the right. This
#  new column contains the slice number, and repeats the series 1 through 16,
#  256 times.

head(baz, 16)

#  [1] 0.94575793 0.05488447 0.81821761 0.66999710
#  [5] 0.40956337 0.63314819 0.12832025 0.14121603
#  [9] 0.04719879 0.25077312 0.96271159 0.67870516
# [13] 0.66355153 0.23132471 0.11800990 0.04486127    
  

Эти первые 16 элементов содержат значения 16 фрагментов в местоположении, указанном в первой строке fooloc . Они идентичны значениям, заданным foo[fooloc[1, 1], fooloc[1, 2], ] .

Полный baz вектор содержит 4096 элементов, соответствующих значениям из 16 слоев, извлеченных в 256 отдельных местоположениях.

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

 qux <- cbind(xyz, baz)
head(qux)

#                       baz
# [1,] 326 264 1 0.94575793
# [2,] 326 264 2 0.05488447
# [3,] 326 264 3 0.81821761
# [4,] 326 264 4 0.66999710
# [5,] 326 264 5 0.40956337
# [6,] 326 264 6 0.63314819
  

Сравнение времени с вашими существующими функциями несправедливо, поскольку pull1 и pull2 возвращает по 32 элемента каждый, и ваш вопрос предполагает, что вам нужно полное 4096. Тем не менее, это все равно быстрее, чем pull1 .

В качестве дополнительного примечания, если вы можете прочитать файл с raster() помощью (из raster пакета), вы можете использовать extract его для извлечения значений из всех слоев, просто указав ваш fooloc в качестве y аргумента.