このページは前リリースの情報です。該当の英語のページはこのリリースで削除されています。

並列木ウェーブレット変換

この例では、信号処理、イメージ処理およびボリューム処理での大きくサンプリングされている DWT に対する並列木複素数離散ウェーブレット変換 (DT-CWT) の利点を示します。並列木 DWT は 2 つの別個の 2 チャネルのフィルター バンクとして実装されます。この例で説明する利点を得るには、2 つのツリーで使用されるスケーリングおよびウェーブレット フィルターを任意に選択することはできません。一方のツリーのローパス (スケーリング) およびハイパス (ウェーブレット) フィルター {h0,h1} は、他方のツリーのローパスおよびハイパス フィルター {g0,g1} によって生成されたスケーリング関数およびウェーブレットの近似ヒルベルト変換である、スケーリング関数およびウェーブレットを生成しなければなりません。したがって、2 つのツリーから形成される複素数値のスケーリング関数およびウェーブレットはほぼ解析的です。

この結果、d 次元データについての冗長性係数が 2d のみである大きくサンプリングされた DWT と比較して、並列木 DWT はシフト変化が少なくなり、方向選択性が強くなります。並列木 DWT の冗長性は、非間引き (定常) DWT の冗長性よりも大幅に少なくなります。

この例では、並列木 DWT の近似シフト不変、ウェーブレットを 2 次元および 3 次元で解析する並列木の選択的配向、およびイメージとボリュームのノイズ除去における並列木の複素数ウェーブレット変換の使用方法について説明します。

並列木 DWT の近接シフト不変

DWT はシフト変化の影響を受けます。つまり、入力信号またはイメージ内の小さなシフトが、DWT 係数のスケール全体における信号/イメージ エネルギーの分布に大きな変化を引き起こします。DT-CWT はほぼシフト不変です。

これをテスト信号で示すには、2 つのシフト離散時間インパルスを 128 サンプルの長さで構築します。一方の信号には 60 番目のサンプルに単位インパルスがありますが、他方の信号には 64 番目のサンプルに単位インパルスがあります。両方の信号に単位エネルギー (2 ノルム) があることは明確です。

kronDelta1 = zeros(128,1);
kronDelta1(60) = 1;
kronDelta2 = zeros(128,1);
kronDelta2(64) = 1;

長さ 14 のウェーブレット フィルターおよびスケーリング フィルターをもつ 2 つの信号の DWT および並列木 DWT をレベル 3 まで下げて取得します。比較のためにレベル 3 の Detail 係数を抽出します。

J = 3;   
dwt1 = dddtree('dwt',kronDelta1,J,'sym7');
dwt2 = dddtree('dwt',kronDelta2,J,'sym7');
dwt1Cfs = dwt1.cfs{J};
dwt2Cfs = dwt2.cfs{J};
      
dt1 = dddtree('cplxdt',kronDelta1,J,'dtf3');
dt2 = dddtree('cplxdt',kronDelta2,J,'dtf3');      
dt1Cfs = dt1.cfs{J}(:,:,1)+1i*dt1.cfs{J}(:,:,2);
dt2Cfs = dt2.cfs{J}(:,:,1)+1i*dt2.cfs{J}(:,:,2);

レベル 3 における 2 つの信号の DWT 係数および DT-CWT 係数の絶対値をプロットし、係数のエネルギー (二乗した 2 ノルム) を計算します。

figure
subplot(1,2,1)
stem(abs(dwt1Cfs),'markerfacecolor',[0 0 1])
title({'DWT';['Squared 2-norm = ' num2str(norm(dwt1Cfs,2)^2,3)]},...
    'fontsize',10)
subplot(1,2,2)
stem(abs(dwt2Cfs),'markerfacecolor',[0 0 1])
title({'DWT';['Squared 2-norm = ' num2str(norm(dwt2Cfs,2)^2,3)]},...
    'fontsize',10)

figure
subplot(1,2,1)
stem(abs(dt1Cfs),'markerfacecolor',[0 0 1])
title({'Dual-tree DWT';['Squared 2-norm = ' num2str(norm(dt1Cfs,2)^2,3)]},...
    'fontsize',10)
subplot(1,2,2)
stem(abs(dwt2Cfs),'markerfacecolor',[0 0 1])
title({'Dual-tree DWT';['Squared 2-norm = ' num2str(norm(dt2Cfs,2)^2,3)]},...
    'fontsize',10)

信号の 4 つのサンプル シフトにより、レベル 3 の DWT ウェーブレット係数のエネルギーに約 6.5% の変化が発生したことに注意してください。ただし、並列木ウェーブレット係数は、エネルギーの変化がたった 0.3% であることを示します。

近似シフト不変の実用性を実際のデータで示すには、心電図 (ECG) 信号を解析します。ECG 信号のサンプリング間隔は 1/180 秒です。データは Percival & Walden[3]の p.125 から採用しています (元々 William Constantine and Per Reinhall, University of Washington によって提供されたデータ)。便宜上、t=0 から始まるデータを採用しました。

load wecg
dt = 1/180;
t = 0:dt:(length(wecg)*dt)-dt;
figure
plot(t,wecg)
xlabel('Seconds')
ylabel('Millivolts')

約 0.7 秒離れている大きい正のピークは、心調律の R 波です。最初に、大きくサンプリングされている DWT を使用する信号を分解し、元の信号とレベル 2 およびレベル 3 のウェーブレット係数をプロットします。R 波が、与えられたサンプリング レートのこれらのスケールで最も著しく分離されているため、レベル 2 およびレベル 3 の係数が選択されました。

J = 6; 
dtDWT1 = dddtree('dwt',wecg,J,'farras');
details = zeros(2048,3);
details(2:4:end,2) = dtDWT1.cfs{2};
details(4:8:end,3) = dtDWT1.cfs{3};
subplot(3,1,1)
stem(t,details(:,2),'Marker','none','ShowBaseline','off')
title('Level 2')
ylabel('mV')
subplot(3,1,2)
stem(t,details(:,3),'Marker','none','ShowBaseline','off')
title('Level 3')
ylabel('mV')
subplot(3,1,3)
plot(t,wecg)
title('Original Signal')
xlabel('Seconds')
ylabel('mV')

並列木変換について上記の解析を繰り返します。この場合、レベル 2 およびレベル 3 での並列木係数の実数部のプロットだけを行います。

dtcplx1 = dddtree('cplxdt',wecg,J,'dtf3');
details = zeros(2048,3);
details(2:4:end,2) = dtcplx1.cfs{2}(:,1,1)+1i*dtcplx1.cfs{2}(:,1,2);
details(4:8:end,3) = dtcplx1.cfs{3}(:,1,1)+1i*dtcplx1.cfs{3}(:,1,2);
subplot(3,1,1)
stem(t,real(details(:,2)),'Marker','none','ShowBaseline','off')
title('Level 2')
ylabel('mV')
subplot(3,1,2)
stem(t,real(details(:,3)),'Marker','none','ShowBaseline','off')
title('Level 3')
ylabel('mV')
subplot(3,1,3)
plot(t,wecg)
title('Original Signal')
xlabel('Seconds')
ylabel('mV')

大きくサンプリングされているウェーブレット変換および並列木ウェーブレット変換の両方によって、同様のスケールに対する ECG の波形の重要な特徴の位置が特定されます。

1 次元信号におけるウェーブレットの重要な応用は、スケールによる分散分析を取得することです。この分散分析は入力信号の循環シフトの影響を受けるべきではないのは当然です。残念ながら、これは大きくサンプリングされている DWT には当てはまりません。これを示すには、ECG 信号を 4 サンプルずつ循環的にシフトして、大きくサンプリングされている DWT でシフトされていない信号とシフトされた信号を解析し、スケール全体のエネルギーの分布を計算します。

wecgShift = circshift(wecg,4);
dtDWT2 = dddtree('dwt',wecgShift,J,'farras');
sigenrgy = norm(wecg,2)^2;
enr1 = cell2mat(cellfun(@(x)(norm(x,2)^2/sigenrgy)*100,dtDWT1.cfs,'uni',0));
enr2 = cell2mat(cellfun(@(x)(norm(x,2)^2/sigenrgy)*100,dtDWT2.cfs,'uni',0));
levels = {'D1';'D2';'D3';'D4';'D5';'D6';'A6'};
enr1 = enr1(:);
enr2 = enr2(:);
table(levels,enr1,enr2,'VariableNames',{'Level','enr1','enr2'})
ans=7×3 table
    Level      enr1      enr2 
    ______    ______    ______

    {'D1'}    4.1994    4.1994
    {'D2'}     8.425     8.425
    {'D3'}    13.381    10.077
    {'D4'}    7.0612    10.031
    {'D5'}    5.4606    5.4436
    {'D6'}    3.1273    3.4584
    {'A6'}    58.345    58.366

元の信号とシフトされた信号の間のレベル 3 および 4 のウェーブレット係数によるエネルギーの変化は約 3% であることに注意してください。次に、複素数並列木ウェーブレット変換を使用してこの解析を繰り返します。

dtcplx2 = dddtree('cplxdt',wecgShift,J,'dtf3');
cfs1 = cellfun(@squeeze,dtcplx1.cfs,'uni',0);
cfs2 = cellfun(@squeeze,dtcplx2.cfs,'uni',0);
cfs1 = cellfun(@(x) complex(x(:,1),x(:,2)),cfs1,'uni',0);
cfs2 = cellfun(@(x) complex(x(:,1),x(:,2)),cfs2,'uni',0);
dtenr1 = cell2mat(cellfun(@(x)(norm(x,2)^2/sigenrgy)*100,cfs1,'uni',0));
dtenr2 = cell2mat(cellfun(@(x)(norm(x,2)^2/sigenrgy)*100,cfs2,'uni',0));
dtenr1 = dtenr1(:);
dtenr2 = dtenr2(:);
table(levels,dtenr1,dtenr2, 'VariableNames',{'Level','dtenr1','dtenr2'})
ans=7×3 table
    Level     dtenr1    dtenr2
    ______    ______    ______

    {'D1'}     4.936     4.936
    {'D2'}    6.6691    6.6691
    {'D3'}    12.682    12.611
    {'D4'}    8.3891    8.4808
    {'D5'}    5.8868    5.8791
    {'D6'}     3.053    3.0415
    {'A6'}    58.384    58.382

並列木変換により、元の信号およびその循環的にシフトされた信号のスケールによる一貫した分散分析が得られます。

画像処理における方向選択性

2 次元における DWT の標準実装では、イメージの列と行の分割可能なフィルター処理を使用します。4 つの消失モーメント (sym4) を持つ Daubechies の最小非対称位相ウェーブレットの LH、HL および HH ウェーブレットを次の図に示します。

figure
J = 5;                      
L = 3*2^(J+1);
N = L/2^J;
Y = zeros(L,3*L);
dt = dddtree2('dwt',Y,J,'sym4');
dt.cfs{J}(N/3,N/2,1) = 1;
dt.cfs{J}(N/2,N/2+N,2) = 1;
dt.cfs{J}(N/2,N/2+2*N,3) = 1;
dwtImage = idddtree2(dt);
imagesc(dwtImage)
axis xy
axis off
title({'Critically Sampled DWT';'2-D separable wavelets (sym4) -- LH, HL, HH'})

LH および HL ウェーブレットにはそれぞれ明確な水平方向および垂直方向があることに注意してください。ただし、右端の HH ウェーブレットは +45 および -45 度の両方の方向を混合して、チェッカーボードのアーティファクトを生成します。この方向の混合は、実数値の分割可能なフィルターの使用によるものです。HH の実数値の分割可能なフィルターには、2 次元周波数平面における高周波数の 4 つの隅すべてに通過帯域があります。

並列木 DWT は、ほぼ解析的である (つまり周波数軸の半分のみをサポートしている) ウェーブレットを使用して方向選択性を実現します。並列木 DWT では、実数部と虚数部の両方に 6 つのサブバンドがあります。6 つの実数部は、2 つのツリーに、列フィルター処理の出力の後に入力イメージの行フィルター処理の出力を追加することで形成されます。6 つの虚数部は、列フィルター処理の出力の後に行フィルター処理の出力を引くことで形成されます。

列と行に適用されるフィルターは、同じフィルターのペア ({h0,h1} または {g0,g1}) にすることも、異なるフィルターのペア ({h0,g1},{g0,h1}) にすることもできます。次のコードは、複素数方向の並列木 DWT の実数部と虚数部に対応する 12 のウェーブレットの方向を示します。

J = 4;
L = 3*2^(J+1);
N = L/2^J;
Y = zeros(2*L,6*L);
wt = dddtree2('cplxdt',Y,J,'dtf3');
wt.cfs{J}(N/2,N/2+0*N,2,2,1) = 1;
wt.cfs{J}(N/2,N/2+1*N,3,1,1) = 1;
wt.cfs{J}(N/2,N/2+2*N,1,2,1) = 1;
wt.cfs{J}(N/2,N/2+3*N,1,1,1) = 1;
wt.cfs{J}(N/2,N/2+4*N,3,2,1) = 1;
wt.cfs{J}(N/2,N/2+5*N,2,1,1) = 1;
wt.cfs{J}(N/2+N,N/2+0*N,2,2,2) = 1;
wt.cfs{J}(N/2+N,N/2+1*N,3,1,2) = 1;
wt.cfs{J}(N/2+N,N/2+2*N,1,2,2) = 1;
wt.cfs{J}(N/2+N,N/2+3*N,1,1,2) = 1;
wt.cfs{J}(N/2+N,N/2+4*N,3,2,2) = 1;
wt.cfs{J}(N/2+N,N/2+5*N,2,1,2) = 1;
waveIm = idddtree2(wt);
imagesc(waveIm)
axis off
title('Complex Dual-Tree 2-D Wavelets')

前の図の先頭行は、実数方向の並列木ウェーブレット変換の 6 つの方向のウェーブレットを示します。2 番目の行は虚数部を示します。実数部と虚数部が共に複素数方向の並列木ウェーブレット変換を形成します。実数部と虚数部の向きは同じ方向です。dddtree2'realdt' オプションと共に使用して、実数部のみを使用する実数方向の並列木 DWT を取得できます。実数方向の並列木ウェーブレット変換を使用すると方向選択性を実現できますが、近似シフト不変などの解析的なウェーブレットを使用するメリットを完全には得られません。

2 次元でのエッジ表現

複素数並列木ウェーブレットの近似解析および選択的配向により、イメージ内のエッジの表現において、2 次元の標準 DWT より優れたパフォーマンスが得られます。これを説明するために、大きくサンプリングされている 2 次元の DWT および 2 次元の複素数方向の並列木変換を使用して、複数方向の直線および曲線の特異点から構成されるエッジをもつテスト イメージを解析します。最初に、線の特異点から構成される八角形のイメージを解析します。

load woctagon
figure
imagesc(woctagon)
colormap gray
title('Original Image')
axis equal
axis off

イメージをレベル 4 まで分解し、レベル 4 の Detail 係数に基づいてイメージの Approximation を再構築します。

dtcplx = dddtree2('cplxdt',woctagon,4,'dtf3');
dtDWT = dddtree2('dwt',woctagon,4,'farras');

dtcplx.cfs{1} = zeros(size(dtcplx.cfs{1}));
dtcplx.cfs{2} = zeros(size(dtcplx.cfs{2}));
dtcplx.cfs{3} = zeros(size(dtcplx.cfs{3}));
dtcplx.cfs{5} = zeros(size(dtcplx.cfs{5}));

dtDWT.cfs{1} = zeros(size(dtDWT.cfs{1}));
dtDWT.cfs{2} = zeros(size(dtDWT.cfs{2}));
dtDWT.cfs{3} = zeros(size(dtDWT.cfs{3}));
dtDWT.cfs{5} = zeros(size(dtDWT.cfs{5}));

dtImage = idddtree2(dtcplx);
dwtImage = idddtree2(dtDWT);
subplot(1,2,1)
imagesc(dtImage)
axis equal
axis off
colormap gray
title('Complex Oriented Dual-Tree')
subplot(1,2,2)
imagesc(dwtImage)
axis equal
axis off
colormap gray
title('DWT')

次に、双曲線エッジをもつ八角形を解析します。双曲線八角形のエッジは曲線の特異点です。

load woctagonHyperbolic
figure
imagesc(woctagonHyperbolic)
colormap gray
title('Octagon with Hyperbolic Edges')
axis equal
axis off

再び、イメージをレベル 4 まで分解し、2 次元の標準 DWT および複素数方向の並列木 DWT の両方についてレベル 4 の Detail 係数に基づいてイメージの Approximation を再構築します。

dtcplx = dddtree2('cplxdt',woctagonHyperbolic,4,'dtf3');
dtDWT = dddtree2('dwt',woctagonHyperbolic,4,'farras');

dtcplx.cfs{1} = zeros(size(dtcplx.cfs{1}));
dtcplx.cfs{2} = zeros(size(dtcplx.cfs{2}));
dtcplx.cfs{3} = zeros(size(dtcplx.cfs{3}));
dtcplx.cfs{5} = zeros(size(dtcplx.cfs{5}));

dtDWT.cfs{1} = zeros(size(dtDWT.cfs{1}));
dtDWT.cfs{2} = zeros(size(dtDWT.cfs{2}));
dtDWT.cfs{3} = zeros(size(dtDWT.cfs{3}));
dtDWT.cfs{5} = zeros(size(dtDWT.cfs{5}));

dtImage = idddtree2(dtcplx);
dwtImage = idddtree2(dtDWT);
subplot(1,2,1)
imagesc(dtImage)
axis equal
axis off
colormap gray
title('DT-CWT')
subplot(1,2,2)
imagesc(dwtImage)
axis equal
axis off
colormap gray
title('DWT')

2 次元の大きくサンプリングされている DWT で明確なリンギング アーティファクトが、両方のイメージの 2 次元の DT-CWT にはないことに注意してください。DT-CWT は、線および曲線の特異点をより忠実に再現します。

イメージのノイズ除去

別個のサブバンド内の異なる方向を分離する機能のため、並列木 DWT はイメージのノイズ除去などの応用において標準の分離可能な DWT より性能が優れていることがよくあります。これを示すために、補助関数 helperCompare2DDenoising を使用します。補助関数はイメージを読み込み、σ=25 であるゼロ平均のホワイト ガウス ノイズを加えます。ユーザー指定のしきい値の範囲の場合、関数は、大きくサンプリングされている DWT、実数方向の並列木 DWT および複素数方向の並列木 DWT のソフトなしきい値処理を使用したノイズ除去を比較します。各しきい値について、平方根平均二乗 (RMS) 誤差およびピーク S/N 比 (PSNR) が表示されます。

numex = 3;
helperCompare2DDenoising(numex,0:2:100,'PlotMetrics');

実数方向および複素数方向の両方の並列木 DWT の性能は、RMS 誤差および PSNR の標準 DWT より優れています。

次に、加法性ノイズの標準偏差と等しい、しきい値 25 のノイズ除去後のイメージを取得します。

numex = 3;
helperCompare2DDenoising(numex,25,'PlotImage');

加法性ノイズの標準偏差と等しいしきい値では、複素数方向の並列木変換は、2 次元の標準 DWT より約 4 dB 高い PSNR を提供します。

3 次元での方向選択性

2 次元の分離可能な DWT で見られるリンギング アーティファクトは、ウェーブレット解析をさらに高い次元に拡張すると悪化します。DT-CWT により、最小限の冗長性で 3 次元の方向選択性を維持できます。3 次元では、並列木変換には 28 のウェーブレット サブバンドがあります。

3 次元の並列木ウェーブレット変換の方向選択性を示すため、3 次元並列木および分離可能な DWT ウェーブレットの両方について 3 次元等値面の例を可視化します。最初に、2 つの並列木サブバンドの実数部と虚数部を別個に可視化します。

helperVisualize3D('Dual-Tree',28,'separate');

helperVisualize3D('Dual-Tree',25,'separate');

等値面プロットの赤い部分は、ウェーブレットのゼロから正の偏位を示し、青い部分は負の偏位を表します。並列木ウェーブレットの実数部と虚数部の空間の方向選択性をはっきりと確認できます。次に、実数部と虚数部のプロットが 1 つの等値面としてまとめてプロットされている並列木サブバンドの 1 つを可視化します。

helperVisualize3D('Dual-Tree',25,'real-imaginary');

前のプロットは、実数部と虚数部が空間内でそれぞれシフトされたものを示します。これは、複素数ウェーブレットの虚数部が実数部の近似ヒルベルト変換であるという事実を反映しています。次に、比較のため、実数の直交ウェーブレットの等値面を 3 次元で可視化します。

helperVisualize3D('DWT',7);

2 次元 DWT で見られる方向の混合は、3 次元ではより顕著になります。2 次元の場合と同様に、3 次元での方向の混合は大幅なリンギングまたはブロッキング アーティファクトを引き起こします。これを示すために、球状ボリュームについて 3 次元の DWT および DT-CWT ウェーブレットの Detail を調べます。球体は 64 × 64 × 64 です。

load sphr
[A,D] = dualtree3(sphr,2,'excludeL1');
A = zeros(size(A));
sphrDTCWT = idualtree3(A,D);
montage(reshape(sphrDTCWT,[64 64 1 64]),'DisplayRange',[])
title('DT-CWT Level 2 Details')

分離可能な DWT に基づいて前のプロットと第 2 レベルの Detail を比較します。

sphrDEC = wavedec3(sphr,2,'sym4','mode','per');
sphrDEC.dec{1} = zeros(size(sphrDEC.dec{1}));
for kk = 2:8
    sphrDEC.dec{kk} = zeros(size(sphrDEC.dec{kk}));
end
sphrrecDWT = waverec3(sphrDEC);
figure
montage(reshape(sphrrecDWT,[64 64 1 64]),'DisplayRange',[])
title('DWT Level 2 Details')

DT-CWT および DWT モンタージュの両方のイメージをズームインすると、DT-CWT と比較して DWT の Detail のブロッキング アーティファクトがどのくらい顕著がわかります。

ボリュームのノイズ除去

2 次元の場合と同様、3 次元の DT-CWT の方向選択性は多くの場合にボリュームのノイズ除去を改善します。

これを示すために、16 スライスで構成される MRI データセットについて考えます。標準偏差 10 のガウス ノイズは元のデータセットに付加されています。ノイズを含むデータセットを表示します。

load MRI3D
montage(reshape(noisyMRI,[128 128 1 16]),'DisplayRange',[])

ノイズ除去前の元の SNR は約 11 dB です。

20*log10(norm(origMRI(:),2)/norm(origMRI(:)-noisyMRI(:),2))
ans = 11.2997

DT-CWT および DWT の両方を使用して、レベル 4 まで下げて MRI データセットのノイズを除去します。同様のウェーブレット フィルター長は両方の場合に使用されます。しきい値の関数として結果の SNR をプロットします。最適な SNR で取得された DT-CWT および DWT の両方についてノイズ除去された結果を表示します。

[imrecDTCWT,imrecDWT] = helperCompare3DDenoising(origMRI,noisyMRI);

figure
montage(reshape(imrecDTCWT,[128 128 1 16]),'DisplayRange',[])
title('DT-CWT Denoised Volume')

figure
montage(reshape(imrecDWT,[128 128 1 16]),'DisplayRange',[])
title('DWT Denoised Volume')

まとめ

並列木 DWT には、大きくサンプリングされている DWT で達成不可能な近接シフト不変および方向選択性の望ましい性質があることを示しました。これらの性質による信号解析のパフォーマンスの向上、イメージとボリューム内でのエッジの表現、およびイメージとボリュームのノイズ除去について説明しました。

参照

[1] Chen, H., and N. G. Kingsbury. "Efficient registration of nonrigid 3-D bodies." IEEE Transactions on Image Processing. Vol. 21, Number 1, January 2012, pp. 262–272.

[2] Kingsbury, N. G. "Complex Wavelets for Shift Invariant Analysis and Filtering of Signals." Journal of Applied and Computational Harmonic Analysis. Vol. 10, Number 3, May 2001, pp. 234–253.

[3] Percival, D. B., and A. T. Walden. Wavelet Methods for Time Series Analysis. Cambridge, UK: Cambridge University Press, 2000.

[4] Selesnick, I., Baraniuk, R. G., and N. G. Kingsbury. "The Dual-Tree Complex Wavelet Transform." IEEE Signal Processing Magazine. Vol. 22, Number 6, November 2005, pp. 123–151.

[5] Selesnick, I. "The Double Density DWT." Wavelets in Signal and Image Analysis: From Theory to Practice (A. A. Petrosian, and F. G. Meyer, eds.). Norwell, MA: Kluwer Academic Publishers, 2001, pp. 39–66.

[6] Selesnick, I. "The Double-Density Dual-Tree Wavelet Transform." IEEE Transactions on Signal Processing. Vol. 52, Number 5, May 2004, pp. 1304–1314.

参考

| | | | |

関連する例

詳細