高速2次元Haar wavelet変換とデータ圧縮

Masahiro MIZUTANI
15 Dec 2011

技術的な基本仮定 N=(2^n)
データ数がN×N個の2次元配列データ signal={s[1,1], s[1,2], ..., s[N,N]} を与えて、高速2次元Haar wavelet変換とその逆Haar wavelet変換を考える。

1-Step 2次元Haar wavelet変換

大きさN=(2^scale)の正方配列  (2^scale)× (2^scale)個の行列標本 データから、1段階のWavelet変換（和平均と平均からのズレ）を係数行列として求める。
係数行列 (N/2)×(N/2)の4ブロックwaveletデータ (Subscript[LL, 1]), (Subscript[LH, 1]),  (Subscript[HL, 1]) (Subscript[HH, 1])からなり、左上の配列ブロック (Subscript[LL, 1])はレベル1のHaar平均行列である。

Haar2DWaveletTransform[signal_List]:=Module[{datasize,scale,hdata,vdata,m,n},
datasize=Length[signal];
(*scale=Log[2,datasize];*)
hdata=ConstantArray[0,{datasize,datasize}];
For[m=0,m≤datasize-1,m++,(* 各行についてWavelet変換 *)
For[n=0,n≤datasize/2-1,n++,
hdata[[m+1,n+1]]=(signal[[m+1,2n+1]]+signal[[m+1,2n+2]])/2;
hdata[[m+1,datasize/2+n+1]]=(signal[[m+1,2n+1]]-signal[[m+1,2n+2]])/2;
(*Print["hdata[",m+1,",",n+1,"]=",hdata[[m+1,n+1]],
"; hdata[",m+1,",",datasize/2+n+1,"]=",hdata[[m+1,datasize/2+n+1]]]*)
]
];
vdata=hdata;
For[n=0,n≤datasize/2-1,n++,(* 次いで各列についてWavelet変換 *)
For[m=0,m≤datasize/2-1,m++,
vdata[[m+1,n+1]]=(hdata[[2m+1,n+1]]+hdata[[2m+2,n+1]])/2;
vdata[[datasize/2+m+1,n+1]]=(hdata[[2m+1,n+1]]-hdata[[2m+2,n+1]])/2;
vdata[[m+1,datasize/2+n+1]]=(hdata[[2m+1,datasize/2+n+1]]+hdata[[2m+2,datasize/2+n+1]])/2;
vdata[[datasize/2+m+1,datasize/2+n+1]]=(hdata[[2m+1,datasize/2+n+1]]-hdata[[2m+2,datasize/2+n+1]])/2;
(*Print["vdata[",m+1,",",n+1,"]=",vdata[[m+1,n+1]],
"; vdata[",datasize/2+m+1,",",n+1,"]=",vdata[[datasize/2+m+1,n+1]],
"; vdata[",m+1,",",datasize/2+n+1,"]=",vdata[[m+1,datasize/2+n+1]],
"; vdata[",datasize/2+m+1,",",datasize/2+n+1,"]=",vdata[[datasize/2+m+1,datasize/2+n+1]]]*)
];
];
vdata
]

t={{9,7,6,2},{5,3,4,4},{8,2,4,0},{6,0,2,2}};

t//MatrixForm

h=Haar2DWaveletTransform[t]

h//MatrixForm

1-Step 2次元Haar wavelet逆変換

Haar2DInverseWaveletTransform[vdata_List]:=Module[{datasize,maxlevel,signal,m,n},
datasize=Length[vdata];
(*maxlevel=Log[2,datasize];*)
signal=ConstantArray[0,{datasize,datasize}];
For[n=0,n≤datasize/2-1,n++,
For[m=0,m≤datasize/2-1,m++,
signal[[2m+1,2n+1]]=vdata[[m+1,n+1]]+vdata[[datasize/2+m+1,n+1]]+vdata[[m+1,datasize/2+n+1]]+vdata[[datasize/2+m+1,datasize/2+n+1]];
signal[[2m+2,2n+1]]=vdata[[m+1,n+1]]-vdata[[datasize/2+m+1,n+1]]+vdata[[m+1,datasize/2+n+1]]-vdata[[datasize/2+m+1,datasize/2+n+1]];
signal[[2m+1,2n+2]]=vdata[[m+1,n+1]]+vdata[[datasize/2+m+1,n+1]]-vdata[[m+1,datasize/2+n+1]]-vdata[[datasize/2+m+1,datasize/2+n+1]];
signal[[2m+2,2n+2]]=vdata[[m+1,n+1]]-vdata[[datasize/2+m+1,n+1]]-vdata[[m+1,datasize/2+n+1]]+vdata[[datasize/2+m+1,datasize/2+n+1]];
]
];
signal
]

Haar2DInverseWaveletTransform[h]

例：画像データ

ImportしたGrayscale画像のデータは0から1までの非負実数である。これをwavelet変換したwavelet係数は負にもなるために、変換した2次元リスト（行列）をそのままImage[]関数を表示しようとしても使っても、うまく表示できない。
下の例では、1回のwavelet変換で求めた係数から定める四区画の行列を画像表示している

SetDirectory[NotebookDirectory[]]

grayconstruct=ImageData[Import["data/grayconst.jpg"]]

tgray=Haar2DWaveletTransform[grayconstruct];

ArrayPlot[tgray]

Image[tgray]

平均画像データを取り出して、元のサイズで表示する

N=512として元のサイズN x N の行列の半分のサイズ (Underoverscript[-, 2, N]) x (Underoverscript[-, 2, N]) の左上行列部分 {1, 256}×{1, 256} には1回Haar wavelet変換によるHaar平均画像データが格納されているので、これを取り出して元のサイズ一片Nで表示する

Image[Take[tgray,{1,256},{1,256}],ImageSize->512]

2次元Haar wavelet変換

1回の2次元Haar変換 Haar2DWaveletTransform[] を繰り返してすべてのWavelet係数を求める

リストデータ lst[m ;; n] = で m番目からn番目までのリストを相当するサイズのリストデータで置き換えるMathematicaの技法を使う

HaarWavelet2DTransform[signal_List]:=Module[{datasize,maxlevel,data,tdata,j},
datasize=Length[signal];
maxlevel=Log[2,datasize];
data=signal;
data=Haar2DWaveletTransform[data];
For[j=1,j<maxlevel,j++,
tdata=Take[data,{1,datasize/2^j},{1,datasize/2^j}];
data[[1;;datasize/2^j,1;;datasize/2^j]]=Haar2DWaveletTransform[tdata];
];
data
]

t={{124,46,47,54,40,43,48,32},{51,45,56,50,40,65,77,40},{32,50,65,44,44,63,65,45},{35,81,55,40,42,60,75,63},{53,69,50,52,65,93,117,104},{45,65,82,74,64,73,75,56},{41,53,68,69,48,40,65,40},{44,57,56,51,44,35,67,58}};

h=HaarWavelet2DTransform[t]

2次元Haar wavelet逆変換

1回の2次元Haar wavelet変換で得たWavelet係数から、1ステップ変換逆関数 Haar2DInverseWaveletTransform[] を繰り返して逆Haar wavelet変換する

HaarInverseWavelet2DTransform[vdata_List]:=Module[{size,scale,wdata,idata,j},
size=Length[vdata];
scale=Log[2,size];
wdata=vdata;
For[j=1,j≤scale,j++,
idata=Take[wdata,{1,2^j},{1,2^j}];
idata=Haar2DInverseWaveletTransform[idata];
wdata[[1;;2^j,1;;2^j]]=idata;
];
wdata
]

HaarInverseWavelet2DTransform[h]

画像データの例

pgrayconstruct=Take[grayconstruct,{256,383},{1,128}];

Length[First[pgrayconstruct]]

Image[pgrayconstruct]

逆変換すると元の画像が表示される

Image[HaarInverseWavelet2DTransform[HaarWavelet2DTransform[pgrayconstruct]]]

データ圧縮

サイズ (2^n)の正方画像圧縮の方法として、指定レベル j のHarr平均データによる圧縮（データ数は (2^n) × (2^n)=(2^(2n))から (2^(n-j)) × (2^(n-j))=(2^(2(n-j)))となって、圧縮率は (2^(-2j))である）と、Wavelet変換係数において微少なwavelet係数を無視し得るとして 0 におきかえ、これを圧縮データとして逆変換して「再生」する方法を考える。

指定した近似レベルの平均化画像

(2^scale)×(2^scale)行列データを1 ステップ変換/逆関数 Haar2DWaveletTransform[] /HaarInverseWavelet2DTransform[] を使って、指定したレベル p までHaar wavelet変換を行ったときに得られたWavelet係数行列の左上にあるサイズ(2^(scale-p))×(2^(scale-p)) 行列はレベルpのHaar wavelet平均である。これを元のサイズで表示する

AveragedHaarWaveletImage[signal, level]
   signal: 元画像の配列データ (前提 (2^n)×(2^n)pixel)
   level: 平均化レベル(1≤level≤n)

AveragedHaarWaveletImage[signal_List,level_Integer]:=Module[{size,scale,data,tdata,j},
size=Length[signal];
scale=Log[2,size];
data=signal;
data=Haar2DWaveletTransform[data];
For[j=1,j<level,j++,
tdata=Take[data,{1,size/2^j},{1,size/2^j}];
data[[1;;size/2^j,1;;size/2^j]]=Haar2DWaveletTransform[tdata];
];
Image[Take[data,{1,size/2^level},{1,size/2^level}],ImageSize->2^scale]
]

AveragedHaarWaveletImage[grayconstruct,3]

微少Wavelet変換係数を無視して0で置き換えて圧縮し、逆変換する

最大振幅に比べて相対的に微少なデータを別な記号に置き換える

ModifiedDataNeglectingRelativeSmallElement[data, thresholdRatio, replacedSymbol]
   [引数] data: 数リスト,   thresholdRatio: 閾値の割合(実数),   replacedSymbol: 置き換え記号
   [戻り値] リストdata内の要素 x の絶対値が、data内の最大絶対値 maxabsの thresholdRaio*maxabs以下 Abs[x] < thresholdRaio*maxabs のとき、そのデータ値を記号  replacedSymbol で置き換える

ModifiedDataNeglectingRelativeSmallElement[data_List,thresholdRatio_Real,replacedSymbol_]:=Module[{maxabs},
maxabs=Max[Abs[data]];
data /. x_ /;Abs[x]<thresholdRatio*maxabs ->replacedSymbol
]

画像データ

grayconstruct=ImageData[Import["data/grayconst.jpg"]];

tgrayconstruct=HaarWavelet2DTransform[grayconstruct];

datanum=Length[Flatten[tgrayconstruct]]

最大絶対値の5% 未満を0にして圧縮

modified05=ModifiedDataNeglectingRelativeSmallElement[tgrayconstruct,0.05,X];

num=Count[Flatten[modified05],X]

圧縮比

Print["Compression Raio for original volume= ", N[(datanum-num)/datanum]*100,"%"]

修正Waveletデータ

modifiedwavelet05=modified05/. X->0;

Image[HaarInverseWavelet2DTransform[modifiedwavelet05]]

最大絶対値の10% 未満を0にして圧縮

modified10=ModifiedDataNeglectingRelativeSmallElement[tgrayconstruct,0.1,X];

num=Count[Flatten[modified10],X]

圧縮比

Print["Compression Raio for original volume= ", N[(datanum-num)/datanum]*100,"%"]

修正Waveletデータ

modifiedwavelet10=modified10/. X->0;

Image[HaarInverseWavelet2DTransform[modifiedwavelet10]]

最大絶対値の15% 未満を0にして圧縮

modified15=ModifiedDataNeglectingRelativeSmallElement[tgrayconstruct,0.15,X];

num=Count[Flatten[modified15],X]

圧縮比

Print["Compression Raio for original volume= ", N[(datanum-num)/datanum]*100,"%"]

修正Waveletデータ

modifiedwavelet15=modified15/. X->0;

Image[HaarInverseWavelet2DTransform[modifiedwavelet15]]

最大絶対値の20% 未満を0にして圧縮

modified20=ModifiedDataNeglectingRelativeSmallElement[tgrayconstruct,0.2,X];

num=Count[Flatten[modified20],X]

圧縮比

Print["Compression Raio for original volume= ", N[(datanum-num)/datanum]*100,"%"]

修正Waveletデータ

modifiedwavelet20=modified20/. X->0;

Image[HaarInverseWavelet2DTransform[modifiedwavelet20]]