Matlab數字圖像處理學習記錄【6】——小波

小波

• 零.前言
• 一.背景知識
• 二. 快速小波變換
• • 2.1 使用小波工具箱的快速小波變換
• 五.圖像進行中的小波
• • 5.1 方向性和邊緣檢測
  • 5.2 圖像平滑或模糊

零.前言

在對數字圖像進行多分辨率觀察和處理時，離散小波變換

DWT

是首選的數學工具。它還能夠幫助我們深入了解圖像的空間域和頻域特性。

一.背景知識

對一個大小為

M×N

的圖像

f(x,y)

，其正向離散變換T(u,v,…)可以用一般的多項式關系表示

T ( u , v , . . . ) = ∑ x , y f ( x , y ) g u , v , . . . ( x , y ) T(u,v,...) = \sum_{x,y}f(x,y)g_{u,v,...}(x,y) T(u,v,...)=x,y∑​f(x,y)gu,v,...​(x,y)

其中，

x和y

為空間變量，而

u,v,...

是變換域變量。若給定T(u,v,…)則有f(x,y)可用的一般的離散反變換得到

f ( x , y ) = ∑ u , v , . . . T ( u , v , . . . ) h u , v , . . . ( x , y ) f(x,y)=\sum_{u,v,...}T(u,v,...)h_{u,v,...}(x,y) f(x,y)=u,v,...∑​T(u,v,...)hu,v,...​(x,y)

其中gu,v和hu,v,…為

正變換核

和

反變換核

。他們決定變換對的

性質

、

計算複雜度

和

主要用途

。

變換系數T(u,v,…)可看做是f關于{hu,v,…}的一系列展開數。

換言之，反變換核對于f的序列展開定義一組展開函數。

DFT

與序列展開表示法完全吻合，這種情況下：

h u , v ( x , y ) = g u , v ∗ ( x , y ) = 1 M N e j 2 π ( u x / M + v y / N ) h_{u,v}(x,y)=g^*_{u,v}(x,y)= \frac{1}{\sqrt{MN}}e^{j2\pi(ux/M+vy/N)} hu,v​(x,y)=gu,v∗​(x,y)=MN

​1​ej2π(ux/M+vy/N)

其中，

j

是虛數機關，

*

是複共轭運算符。該變換域的v和u分别表示

水準頻度

和

垂直頻度

。變換核是可分的：

h u , v ( x , y ) = h u ( x ) h v ( y ) h_{u,v}(x,y)=h_u(x)h_v(y) hu,v​(x,y)=hu​(x)hv​(y)

h u ( x ) = 1 M e j 2 π u x / M 和 h v ( y ) = 1 N e j 2 π v y / N h_u(x) = \frac{1}{\sqrt{M}}e^{j2\pi ux/M}和h_v(y) = \frac{1}{\sqrt{N}}e^{j2\pi vy/N} hu​(x)=M

​1​ej2πux/M和hv​(y)=N

​1​ej2πvy/N

且該兩個分量是正交的，因為：

⟨ h r , h s ⟩ = δ r s = { 1 , r = s 0 , 其 他 \left \langle h_r, h_s\right \rangle = \delta_{rs}= \begin{cases} 1, & r=s\\ 0, & 其他 \end{cases} ⟨hr​,hs​⟩=δrs​={1,0,​r=s其他​

變換核的可分性簡化了二維變換的計算，因為這樣我們就可以對行分開進行一維變換來實作二維變換。正交性則導緻了正反變換核之間的複共轭關系。

與離散傅裡葉變換不同,離散傅裡葉變換完全可以通過關于一對(前面已給出的）變換核的兩個簡單方程來定義,而術語“離散小波變換”指的是這樣一類變換：

不僅其中使用的變換核不同(例如，所用的展開函數)，而且這些函數（例如，不管它們構成正交基還是雙正交基）的基本特性和它們應用的方法（例如，計算多少不同的分辨率)也不同。因為DWT包含各種獨特但相關的變換，是以我們不能寫出一個能完全描述它們的公式。相反,我們利用變換核對或定義該核對的一組參數來表征每個DWT。各種變換都與這樣的事實有關，即變換的展開函數是變化頻率和持續時間受限的“小波”。

傅裡葉展開後，是由無窮個在持續時間無限的正弦波

DWT展開後，是持續時間有限變化的

小波

在本章剩餘的部分，将介紹幾種小波核。每一個小波核都有如下的基本特性：

1. 可分離性、尺度可變性和平移性。核可用三個可分的二維小波表示

   ψ H ( x , y ) = ψ ( x ) φ ( y ) ψ V ( x , y ) = φ ( x ) ψ ( y ) ψ D ( x , y ) = ψ ( x ) ψ ( y ) \psi^H(x,y)=\psi(x)\varphi(y) \\ \psi^V(x,y)=\varphi(x)\psi(y)\\ \psi^D(x,y)=\psi(x)\psi(y) ψH(x,y)=ψ(x)φ(y)ψV(x,y)=φ(x)ψ(y)ψD(x,y)=ψ(x)ψ(y)

   這三者從上到下依次被稱為：水準、垂直和對角小波，且二維可分的尺度函數是

   φ ( x , y ) = φ ( x ) φ ( y ) \varphi(x,y)= \varphi(x)\varphi(y) φ(x,y)=φ(x)φ(y)

   每個二維函數是兩個一維實平方可積的尺度和小波函數的乘積：

   φ j , k ( x ) = 2 j / 2 φ ( 2 j x − k ) ψ j , k ( x ) = 2 j / 2 ψ ( 2 j x − k ) \varphi_{j,k}(x) = 2^{j/2}\varphi(2^jx-k) \\ \psi_{j,k}(x)=2^{j/2}\psi(2^jx-k) φj,k​(x)=2j/2φ(2jx−k)ψj,k​(x)=2j/2ψ(2jx−k)

   平移參數

   k

   決定一位函數沿x軸的位置，尺度

   j

   決定寬度。2^{j/2}控制他們的高度和振幅。
2. 多分辨率的一緻性。剛才介紹的一維尺度函數滿足多分辨率分析的如下需求：

• φj,k與其整數平移正交
• 在低尺度下可表示為一系列φj,k展開的一組函數，包含在可以以更高尺度表示的那些函數中。
• 唯一可以任意尺度表示的函數是f(x)=0

• 當j→∞時，可以用任意精度來表示任意函數。

  當這些條件滿足是，就存在一個伴随小波ψj,k及其整數平移和二進制尺度，其變化範圍是在鄰接尺度上可以表示的任意兩組φj,k之間的差。

1. 正交性。 展開函數（如{φj,k(x)}）對于一組一維可測的、平方可積函數形成一個正交基或雙正交基。之是以稱為基，是因為對于每一個可描述函數必須有惟一一組展開系數。正如在介紹傅裡葉核時所說明的那樣,對于實數來說正交核為gu,v,… = hu,v,…。對于雙正交情況

   ⟨ h r , h s ⟩ = δ r s = { 1 , r = s 0 , 其 他 \left \langle h_r, h_s\right \rangle = \delta_{rs}= \begin{cases} 1, & r=s\\ 0, & 其他 \end{cases} ⟨hr​,hs​⟩=δrs​={1,0,​r=s其他​

   并且g成為h的對偶。對于利用尺度和小波函數φj,k(x)和ψj,k(x)的雙正交小波變換。對偶分别定義為 φ ~ j , k ( x ) 和 ψ ~ j , k ( x ) \tilde{\varphi}_{j,k}(x)和\tilde{\psi}_{j,k}(x) φ~​j,k​(x)和ψ~​j,k​(x)

說了那麼多小波的性質，個人了解為：通過一種固定的信号(

小波核

)，且該信号在時域上均值為0的，去“比對”給定信号。通過二分法縮放其持續時間，來修改其頻率。而又通過平移來修改小波在時域上“波動的位置”。這樣能夠更好的提取某個信号的某時段的頻率資訊，而不是單純的像傅裡葉那樣得到一個頻譜。感覺就像是将在變換函數中，加入了可調節視窗大小的特殊的傅裡葉變換。

二. 快速小波變換

上面性質的一個重要結果是φ(x)和ψ(x)可以用它們自身的雙分辨率副本的線性組合來表達。這樣序列展開後：

φ ( x ) = ∑ n h φ ( n ) 2 φ ( 2 x − n ) ψ ( x ) = ∑ n h ψ ( n ) 2 φ ( 2 x − n ) \varphi(x) = \sum_nh_\varphi(n)\sqrt{2}\varphi(2x-n) \\ \psi(x) = \sum_nh_\psi(n)\sqrt{2}\varphi(2x-n) φ(x)=n∑​hφ​(n)2

​φ(2x−n)ψ(x)=n∑​hψ​(n)2

​φ(2x−n)

其中，hφ和hψ展開系數分别稱為尺度和小波向量。它們是快速小波變換

FWT

濾波器的系數。DWT的疊代計算方法如圖所示：

在途中，輸出Wφ(j,m,n)和{Wψi(j,m,n)，i=H,V,D}是在遲鈍j處的DWT系數。方框中包含了時間反轉尺度和小波向量。單元hφ(-n)和hψ(-m)分别是低通和高通分解濾波器。最後最後，包括2和向下箭頭的方框表示下取樣，即從點的序列中每隔一個點來提取一個點。準确地說，用于計算圖7.2中的WψH(i, m, n)的濾波和下取樣操作是：

W ψ H ( j , m , n ) = h ψ ( − m ) ∗ [ h φ ( − n ) ∗ W φ ( j + 1 , m , n ) ∣ n = 2 k , k ≥ 0 ] ∣ m = 2 k , k ≥ 0 W_\psi^H(j,m,n)=h_\psi(-m)*[h_\varphi(-n)*W_\varphi(j+1,m,n)|_{n=2k,k≥0}]|_{m=2k, k≥0} WψH​(j,m,n)=hψ​(−m)∗[hφ​(−n)∗Wφ​(j+1,m,n)∣n=2k,k≥0​]∣m=2k,k≥0​

其中，

*

表示卷積。在非負處計算卷積時，偶數索引等價于濾波和用2下取樣。

2.1 使用小波工具箱的快速小波變換

小波工具箱對于各種快速小波變換提供了分解濾波器。與特殊變換有關的濾波器可以通過

wfilters

來通路：

[Lo_D, Hi_D, Lo_R, Hi_R] = wfilters(wname)

其中，輸入參數

wname

決定傳回的濾波器系數。輸出四個行向量，他們分别是低通分解、高通分解、低通重構、高通重構濾波器。

頻繁地耦合的濾波器可以用如下語句交替檢索：

[F1, F2] = wfilters(wname, type)

将type設為

d、r、l、h

可以分别或得一堆分解、重構、低通和高通濾波器。被分解或低通濾波器在F1裡傳回，重構或高通放在F2中。

上表列出了包含在小波工具箱中的FWT濾波器。它們的性質以及其他與尺度和小波函數相關的資訊可在數字濾波和多分辨率分析方面的文獻中找到。一些更重要的特性由小波工具箱的函數waveinfo和wavefun提供。例如，為了在MATLAB的指令視窗列印小波族wfamily，可輸入

waveinfo(wfamily)

标準正交變換的尺度或小波函數的函數的數字近似隻，可通過：

[phi, psi, xval] = wavefun(wname, iter)

獲得，傳回近似向量φ和ψ并且計算向量xval。參數

iter

控制疊代次數，決定近似值的精度。對于雙正交變換，合适的文法為：

[phi1, psi1, phi2, psi2, xval] = wavefun(wname， iter)

phi1

和

psi1

是分解函數，

phi2

和

psi2

是重構函數

通過以下代碼我們能夠查詢一個小波函數的圖形：

[phi, psi, xval] = wavefun('haar', 10);
xaxis = zeros(size(xval));
subplot(1, 2, 1);
plot(xval, phi, 'k', xval, xaxis, '--k');
axis([0 1 -1.5 1.5]); axis square;
subplot(1, 2, 2);
plot(xval, psi, 'k', xval, xaxis, '--k');
axis([0 1 -1.5 1.5]); axis square;
           

然後可以用兩種方式對函數進行小波變換：

[C S] = wavedec2(X, N, Lo_D, Hi_D)

和

[C S] = wavedec2(X, N, wname)

其中

X

是二維矩陣，

N

是被計算的尺度，後者要麼給出分解後的濾波器，或者是小波的名稱。

C

是由行向量C——小波變換系數

矩陣

S

是定義了C中資料排列的記錄矩陣S

後面就是一些小波變換代碼層面的數學應用，由于時間太碎片化，并沒有太多時間去認真分析這些，是以将該處内容留存至以後再學。

五.圖像進行中的小波

這節則是對圖像處理的小波的應用。基于小波的圖像處理辦法過程因該是：

1. 計算一幅圖像的而為小波變換
2. 修改變換系數
3. 計算反變換

5.1 方向性和邊緣檢測

f = imread('./example.tif');
subplot(2, 2, 1);
imshow(f);
subplot(2, 2, 2);
[c, s] = wavefast(f, 1, 'sym4');
wave2gray(c, s, -6);
subplot(2, 2, 3);
[nc, y] = wavecut('a', c, s);
wave2gray(nc, s, -6);
edge = abs(waveback(nc, s, 'sym4'));
subplot(2, 2, 4);
imshow(mat2gray(edge));
           

關于’sym4’小波，圖1的單尺度小波變換的水準、垂直、對角線的方向性在圖2中可清楚地看到。注意，原圖像的水準邊緣出現在圖2的右上象限的水準細節系數中。對于圖像的垂直邊緣,可以在左下象限的垂直細節系數中類似地确定。要将這些資訊合并成為一幅邊緣圖像，我們可以簡單地把生成的變換的近似系數設為零,計算它的反變換,再對其取絕對值。修改過的變換和得到的邊緣圖像分别顯示在圖3和圖4中。類似的過程可用于隔離垂直邊緣或水準邊緣。

5.2 圖像平滑或模糊

通過

wavezero

函數 将某一級小波變換的值置0.

function [nc, g8] = wavezero(c, s, l, wname)
%WAVEZERO Zeroes wavelet transform detail coefficients. 
%   [NC, G8] = WAVEZERO(C, S, L, WNAME) zeroes the level L detail
%   coefficients in wavelet decomposition structure [C, S] and
%   computes the resulting inverse transform with respect to WNAME
%   wavelets.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.4 $  $Date: 2003/10/13 01:31:35 $

[nc, foo] = wavecut('h', c, s, l);
[nc, foo] = wavecut('v', nc, s, l);
[nc, foo] = wavecut('d', nc, s, l);
i = waveback(nc, s, wname);
g8 = im2uint8(mat2gray(i));
figure; imshow(g8);
           

然後通過不斷調用

wavezero

将圖檔更平滑：

f = imread('./example.tif')
[c, s] = wavefast(f, 4, 'sym4');
wave2gray(c, s, 20);
[c, g8] = wavezero(c, s, 1, 'sym4');
[c, g8] = wavezero(c, s, 2, 'sym4');
[c, g8] = wavezero(c, s, 3, 'sym4');
[c, g8] = wavezero(c, s, 4, 'sym4');
           

然後發現級數越高模糊程度越高。是以，小波函數域中的尺度和傅裡葉域中的頻率之間的有着某種聯系。