2014-11-28

任意次元〜私のための離散ウェーブレット変換

昨日の記事で、ウェーブレット変換をべたべたになぞって、２次元をやってみた
その記事では、変換後の値の定数倍とかは気にせず、各セルの大小比がとれていればよい、という、粗雑なものだが、それでも、２次元にする、というのは、１次元のウェーブレット処理を、縦と横とに順繰りにかけるだけであることも示した
さらに、高次元にするのも同じやり口であることがわかっている
この記事では、任意次元で処理するための関数を作ってみる
１次元のベクトルを２次元の行列に上げ、さらに任意次元のアレイに上げたい
そのために、こちらの記事では、アレイの転置を、そしてこちらの記事では、Rでアレイにapply()関数をかけつつ、アレイのセルアドレスを動かさないための工夫をした
アレイの処理のための関数２つ

my.t.array <- function(A,v){
	A.v <- c(A)
	p <- list()
	as <- dim(A)
	for(i in 1:length(as)){
		p[[i]] <- 1:as[i]
	}
	ad <- as.matrix(expand.grid(p))

	as.2 <- as[v]
	val2 <- apply(v <- (t(ad[,v])-1) * c(1,cumprod(as.2)[-length(as.2)]),2,sum)+1
	array(A.v[order(val2)],as.2)
}
my.apply.fix.dim <- function(A,k,f, ...){
	as <- dim(A)
	d <- length(as)
	no.k <- (1:d)[-k]
	post.apply <- apply(A,k,f,...)
	tmp <-c(no.k,k)
	my.t.array(array(post.apply,as[tmp]),order(tmp))
}

それを使って任意次元のウェーブレット変換

my.dwt.2d <- function(x,ord=c(1,2),fs){
# 1次元データにフィルタfをかけて、xと同じ長さのベクトルを返す関数
	tmp.fx <- function(x,f){
		x. <- c(rep(0,length(f)-1),x)
		y <- rep(0,length(x))
		for(i in 1:length(x)){
			y[i] <- sum(x.[i:(i+length(f)-1)] * f)
		}
		y
	}
	n <- length(x[,1])
	X <- x
# 上で作った関数を、行についてfs[[1]]のフィルタをかけ
# ついで、列についてfs[[2]]のフィルタをかける
	for(i in 1:length(ord)){
		for(j in 1:n){
			if(ord[i] ==1){
				X[j,] <- tmp.fx(X[j,],fs[[i]])
			}else if(ord[i]==2){
				X[,j] <- tmp.fx(X[,j],fs[[i]])
			}
		}
	}
	tmp <- 1:(n/2)
# 1/4個のセルを取り出す
	X[(tmp*2),(tmp*2)]
}
RC.out.f1.f1 <- my.dwt.2d(xbox,c(1,2),list(f1,f1))
RC.out.f1.f2 <- my.dwt.2d(xbox,c(1,2),list(f1,f2))
RC.out.f2.f1 <- my.dwt.2d(xbox,c(1,2),list(f2,f1))
RC.out.f2.f2 <- my.dwt.2d(xbox,c(1,2),list(f2,f2))
# waveslimパッケージの関数を使う
xbox.out <- dwt.2d(xbox,"haar",1)
# ４つずつできた出力行列は同じ色彩(値は定数倍の違いだけ)
par(mfcol=c(2,2))
image(RC.out.f2.f1)
image(RC.out.f1.f2)
image(RC.out.f1.f1)
image(RC.out.f2.f2)
par(mfcol=c(1,1))
dev.new()

par(mfcol=c(2,2))
image(xbox.out[[1]])
image(xbox.out[[2]])
image(xbox.out[[3]])
image(xbox.out[[4]])
par(mfcol=c(1,1))
# xはn次元立方アレイ
# ordはウェーブレット処理する軸の順番
# fsはリストで、離散ウェーブレットの波関数。軸ごとに指定
my.dwt.nd <- function(x,ord=1:length(dim(x)),fs){
# ベクトルのウェーブレット変換のために、fs[[i]]をかぶせる内部関数
	tmp.fx <- function(x,f){
		x. <- c(rep(0,length(f)-1),x)
		y <- rep(0,length(x))
		for(i in 1:length(x)){
			y[i] <- sum(x.[i:(i+length(f)-1)] * f)
		}
		y
	}
# アレイの軸長のベクトルと次元
	as <- dim(x)
	d <- length(dim(x))
	n <- length(x)/(2^d)
	X <- x
# 軸ごとに関数をかぶせて、中身を置き換える
	for(i in 1:length(ord)){
# この関数が、アレイの軸iごとにウェーブレット関数をかぶせつつ、アレイの形をとどめる処理
		X <- my.apply.fix.dim(X,ord[-i],tmp.fx,fs[[i]])
	}
# 2^dセルのうち、１つが取り出される（全軸のセル番号が偶数のものを取り出す
	address <- which((X==0 | X!=0),arr.ind=TRUE)
	s <- apply(address %% 2,1,sum)
	array(X[which(s==0)],dim(x)/2)
}

使ってみる
まずは、２次元で試す

library(waveslim)
k1 <- 4
k2 <- 2
x <- array(1:(2^(k1*k2)),c(2^k1,2^k1))

fs <- list(c(1,1),c(1,1))
# 任意次元関数の出力と、2次元用に作った関数（動作確認済み）との結果の一致を確認
my.out2.nd <- my.dwt.nd(x,fs=fs)
my.out2.2d <- my.dwt.2d(x,fs=fs)
out.dwt2 <- dwt.2d(x,"haar",1)
plot(c(my.out2.nd),c(my.out2.2d))
# waveslimパッケージの２次元用関数の出力との定数倍関係を確認
plot(c(my.out2.nd),c(out.dwt2$LL1))

３次元

k1 <- 3
k2 <- 3
x <- array(rnorm(2^(k1*k2)),c(2^k1,2^k1,2^k1))
fs <- list(c(1,1),c(1,1),c(1,1))
my.out3.nd <- my.dwt.nd(x,fs=fs)
out.dwt3 <- dwt.3d(x,"haar",1)
# c(1,1)はLにc(1,-1)はHに対応するので、いくつかwaveslimパッケージのdwt.3d()の出力と比較して定数倍関係にあることを確認
plot(c(my.out3.nd),c(out.dwt3$LLL1))
fs <- list(c(1,1),c(1,-1),c(1,-1))
my.out3.nd <- my.dwt.nd(x,fs=fs)
out.dwt3 <- dwt.3d(x,"haar",1)

plot(c(my.out3.nd),c(out.dwt3$LHH1))