多次元アレイの畳み込みフィルタ

• こちらで四元数離散フーリエ変換をやっている
• その中で多次元アレイで畳み込みも実施したい
• やってみる
• 多次元アレイとフィルタアレイを与えると、畳み込んでほしい
• 多次元アレイをフィルタの大きさに応じて値が0であるようなセルの縁取りをすることで拡大して、フィルタのすべてのセルについて、その中心からのずれ番地を考慮しながら元の多次元フィルタの値の重みづけ総和を出す

# アレイの番地を出す補助関数
my.array.address <- function(a){
	d <- dim(a)
	L <- list()
	L[[1]] <- 1:d[1]
	for(i in 2:length(d)){
		L[[i]] <- 1:d[i]
	}
	as.matrix(expand.grid(L))
}

a <- array(0,c(2,3,4))
my.array.address(a)
# アレイアドレスを任意の軸別進数にての順序数に変える補助関数
my.array.loc <- function(ad,d){
	d. <- c(1,cumprod(d)[1:(length(d)-1)])
	apply((t(ad)-1) * d.,2,sum)+1
}
# アレイに縁取り幅を指定して大型化する補助関数
# d1 は各次元で1より前に付け加える行数
# d2 は最終行より後に付け加える行数
# d1,d2 < 0 ならば、縮める作業
my.array.expansion <- function(a,d1,d2=d1){
	D <- dim(a)
	ad <- my.array.address(a)
	ad.new <- t(t(ad) + d1)
	D.new <- D + d1 + d2
	ret <- array(0,D.new)
	tmp <- ad.new > 0
	tmp2 <-t(t(ad.new) - D.new) <=0
	s <- which(apply(tmp*tmp2,1,prod)==1)
	loc.new <- my.array.loc(ad.new[s,],D.new)
	ret[loc.new] <- a[s]
	ret
}
d <- c(2,3,4)
a <- array(1:prod(d),d)
d1 <- rep(2,length(d))
d2 <- rep(3,length(d))
d1 <- c(1,1,1)
d2 <- c(2,4,2)
my.array.expansion(a,d1,d2)

# フィルタリング関数
# aは任意次元アレイ
# fはaと同じ軸数を持つフィルタアレイ
# ctrはフィルタの中心アドレス(その値のみが１で残りのセルは全部０のフィルタfの場合、返り値は入力aそのものとなる)
my.array.filter <- function(a,f,ctr=(dim(f)+1)/2){
	d.a <- dim(a)
	d.f <- dim(f)
	diff.d1 <- ctr-1
	diff.d2 <- d.f-ctr
	ad.f <- my.array.address(f)
	ad.f. <- ad.f - ctr
	ad.a <- my.array.address(a)
	D.new <- d.a + diff.d1 + diff.d2
	#a.big <- my.array.expansion(a,diff.d1,diff.d2)
	a.big <- array(0,D.new)
	max.loc <- prod(D.new)
	for(i in 1:length(ad.f.[,1])){
		tmp.ad <- t(t(ad.a) + (-1)*ad.f.[i,]+diff.d1)
		tmp.v <- a * f[i]
		tmp.loc <- my.array.loc(tmp.ad,D.new)
		s <- which(tmp.loc>0 & tmp.loc<max.loc)
		a.big[tmp.loc[s]] <- a.big[tmp.loc[s]] + tmp.v[s]
	}
	my.array.expansion(a.big,-diff.d1,-diff.d2)
}
d <- c(2,3,4)
a <- array(1:prod(d),d)
#f <- array(1,rep(3,length(d)))
f <- array(0,rep(3,length(d)))
#f <- f - 1
#f[2,2] <- 2
f[2,2,2] <- 1
f[1,1,1] <- 1
a. <- my.array.filter(a,f)
a. - a

> a. - a
, , 1

     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0

, , 2

     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    1    3

, , 3

     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    7    9

, , 4

     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0   13   15