数字图像处理MATLAB学习笔记(六)
Wavelets
1 Background
这里我们简单回顾一下傅里叶变换,并对一些概念做一下定义:
设正向DFT为
T
(
u
,
v
,
…
)
{T(u,v,\dots)}
T(u,v,…),根据一般关系DFT和IDFT关系式可以表示为:
T
(
u
,
v
,
…
)
=
∑
x
,
y
f
(
x
,
y
)
g
u
,
v
,
…
(
x
,
y
)
f
(
x
,
y
)
=
∑
x
,
y
T
(
u
,
v
,
…
)
h
u
,
v
,
…
(
x
,
y
)
T(u,v,\dots)=\sum_{x,y}f(x,y)g_{u,v,\dots}(x,y)\\ f(x,y)=\sum_{x,y}T(u,v,\dots)h_{u,v,\dots}(x,y)
T(u,v,…)=x,y∑f(x,y)gu,v,…(x,y)f(x,y)=x,y∑T(u,v,…)hu,v,…(x,y)
其中
g
u
,
v
,
…
{g_{u,v,\dots}}
gu,v,…和
h
u
,
v
,
…
{h_{u,v,\dots}}
hu,v,…分别表示为正变换核(Lowpass)和反变换核(Highpass)。我们用级数展开可以表达关系为:
h
u
,
v
(
x
,
y
)
=
g
u
,
v
∗
(
x
,
y
)
=
1
M
N
exp
(
j
2
π
(
u
x
/
M
+
v
y
/
N
)
)
h_{u,v}(x,y)=g_{u,v}^*(x,y)=\frac{1}{\sqrt{MN}}\exp(j2\pi(ux/M+vy/N))
hu,v(x,y)=gu,v∗(x,y)=MN1exp(j2π(ux/M+vy/N))
其中,
j
=
−
1
{j=\sqrt{-1}}
j=−1,
∗
{*}
∗表示为复共轭算子,
u
=
0
,
1
,
2
,
…
,
M
−
1
{u=0,1,2,\dots,M-1}
u=0,1,2,…,M−1且
v
=
0
,
1
,
2
,
…
,
N
−
1
{v=0,1,2,\dots,N-1}
v=0,1,2,…,N−1。u和v分别表示水平和垂直频率,变换核可分为
h
u
(
x
,
y
)
=
1
M
exp
(
j
2
π
x
/
M
)
h
v
(
x
,
y
)
=
1
N
exp
(
j
2
π
y
/
N
)
且,二者是正交的,正交关系为:
<
h
r
,
h
s
>
=
δ
r
s
=
{
1
r
=
s
0
otherwise
h_u(x,y)=\frac{1}{\sqrt{M}}\exp(j2\pi x/M)\\ h_v(x,y)=\frac{1}{\sqrt{N}}\exp(j2\pi y/N)\\ \text{且,二者是正交的,正交关系为:}\\ \big<h_r,h_s\big>=\delta_{rs}=\begin{cases}1&r=s\\0&\text{otherwise}\end{cases}
hu(x,y)=M1exp(j2πx/M)hv(x,y)=N1exp(j2πy/N)且,二者是正交的,正交关系为:⟨hr,hs⟩=δrs={10r=sotherwise
DWT(Discrete Wavelets Transform)离散小波变换指的是基本小波的尺度和平移进行离散化。
针对小波核,具有的一般特性为:
- 可分离性、可伸缩性、可平移性:
变换核可用三个可分的二维小波来表示:
ψ
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)
其中,
ψ
H
(
x
,
y
)
{\psi^H(x,y)}
ψH(x,y)、
ψ
V
(
x
,
y
)
{\psi^V(x,y)}
ψV(x,y)和
ψ
D
(
x
,
y
)
{\psi^D(x,y)}
ψD(x,y)分别为水平Horizontal、垂直Vertical、对角Diagonal小波。且是二维可分的尺度函数为
φ
(
x
,
y
)
=
φ
(
x
)
φ
(
y
)
{\varphi(x,y)=\varphi(x)\varphi(y)}
φ(x,y)=φ(x)φ(y)。
每个二维函数都是两个一维是函数的、平方可积的尺度和小波函数的乘积: φ j , k = 2 j / 2 φ ( 2 j x − k ) {\varphi_{j,k}=2^{j/2}\varphi(2^jx-k)} φj,k=2j/2φ(2jx−k)和 ψ j , k = 2 j / 2 ψ ( 2 j x − k ) {\psi_{j,k}=2^{j/2}\psi(2^jx-k)} ψj,k=2j/2ψ(2jx−k)。平移参数 k k k决定这些一维函数沿x轴的位置,尺度j决定他们的宽度并且按照 2 j / 2 2^{j/2} 2j/2控制高度或振幅。
其中,与展开函数相关联的是母小波的二进制尺度和证书平移 ψ ( x ) = ψ 0 , 0 ( x ) {\psi(x)=\psi_{0,0}(x)} ψ(x)=ψ0,0(x),并且尺度函数 φ ( x ) = φ 0 , 0 ( x ) {\varphi(x)=\varphi_{0,0}(x)} φ(x)=φ0,0(x)。
- Multiresolution Compatibility多分辨率的兼容性:一维尺度函数满足下要求即可用于多分辨率分析:
- φ j , k {\varphi_{j,k}} φj,k对整数平移是正交的。
- 在低尺度或低分辨率下,以一系列 φ j , k {\varphi_{j,k}} φj,k的展开函数来描述的一组函数包含在可以用高尺度描述的函数中。
- 唯一可在每个尺度上描述的函数是 f ( x ) = 0 f(x)=0 f(x)=0。
- 在 j → ∞ j\to\infty j→∞,任何函数都可以以任意精度来描述。
- Orthogonality正交性
小波变换和Fourier Transform之间的关系:
原始信号中,时域的分辨率高,频域分辨率为零,意味着我们可以在时域中分辨出非常小的特征,而在频域中没有特征。
傅里叶变换在频域具有高分辨率,在时域具有零分辨率。
短时傅里叶变换相对于上面的两种方法,它在时域和频域具有中等的分辨率。
小波变换有如下的特点:
-
对于小频率值,频域分辨率高,时域分辨率低
-
对于大频率值,频域分辨率低,时域分辨率高
小波变换在频域分辨率和时域分辨率两者之间做了权衡:在与时间相关的特征上上,它在时域具有高分辨率;而在与频率相关的特征上,它在频域具有高分辨率。
2 The Fast Wavelets Transform(FWT)
根据小波核,我们有个重要推论: φ ( x ) \varphi(x) φ(x)和 ψ ( x ) \psi(x) ψ(x)可以用双分辨率副本的线性组合来表达。
通过级数展开,我们得到:
φ
(
x
)
=
∑
n
h
φ
(
n
)
2
φ
(
2
x
−
n
)
ψ
(
x
)
=
∑
n
h
ψ
(
n
)
2
φ
(
2
x
−
n
)
\varphi(x)=\sum_nh_\varphi(n)\sqrt2\varphi(2x-n)\\ \psi(x)=\sum_nh_\psi(n)\sqrt2\varphi(2x-n)
φ(x)=n∑hφ(n)2φ(2x−n)ψ(x)=n∑hψ(n)2φ(2x−n)
其中,
h
φ
h_{\varphi}
hφ和
h
ψ
h_\psi
hψ分别是尺度sacling和小波矢量wavelet vectors。
关于FWT的迭代计算如图:
![image-20211201175312011](/Users/wanghe/Library/Application%20Support/typora-user-images/image-20211201175312011.png)
用公式表达如下:
W
ψ
D
(
j
,
m
,
n
)
=
h
ψ
(
−
m
)
⊙
[
h
ψ
(
−
n
)
⊙
W
φ
(
j
+
1
,
m
,
n
)
∣
n
=
2
k
,
k
≥
0
]
∣
m
=
2
k
,
k
≥
0
W
ψ
V
(
j
,
m
,
n
)
=
h
φ
(
−
m
)
⊙
[
h
ψ
(
−
n
)
⊙
W
φ
(
j
+
1
,
m
,
n
)
∣
n
=
2
k
,
k
≥
0
]
∣
m
=
2
k
,
k
≥
0
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
φ
(
j
,
m
,
n
)
=
h
φ
(
−
m
)
⊙
[
h
φ
(
−
n
)
⊙
W
φ
(
j
+
1
,
m
,
n
)
∣
n
=
2
k
,
k
≥
0
]
∣
m
=
2
k
,
k
≥
0
W_\psi^D(j,m,n)=h_\psi(-m)\odot[h_\psi(-n)\odot W_\varphi(j+1,m,n)|_{n=2k,k\ge0}]|_{m=2k,k\ge0}\\ W_\psi^V(j,m,n)=h_\varphi(-m)\odot[h_\psi(-n)\odot W_\varphi(j+1,m,n)|_{n=2k,k\ge0}]|_{m=2k,k\ge0}\\ W_\psi^H(j,m,n)=h_\psi(-m)\odot[h_\varphi(-n)\odot W_\varphi(j+1,m,n)|_{n=2k,k\ge0}]|_{m=2k,k\ge0}\\ W_\varphi(j,m,n)=h_\varphi(-m)\odot[h_\varphi(-n)\odot W_\varphi(j+1,m,n)|_{n=2k,k\ge0}]|_{m=2k,k\ge0}
WψD(j,m,n)=hψ(−m)⊙[hψ(−n)⊙Wφ(j+1,m,n)∣n=2k,k≥0]∣m=2k,k≥0WψV(j,m,n)=hφ(−m)⊙[hψ(−n)⊙Wφ(j+1,m,n)∣n=2k,k≥0]∣m=2k,k≥0WψH(j,m,n)=hψ(−m)⊙[hφ(−n)⊙Wφ(j+1,m,n)∣n=2k,k≥0]∣m=2k,k≥0Wφ(j,m,n)=hφ(−m)⊙[hφ(−n)⊙Wφ(j+1,m,n)∣n=2k,k≥0]∣m=2k,k≥0
其中,
h
φ
h_\varphi
hφ和
h
ψ
h_\psi
hψ分别是lowpass 和 highpass decomposition filters。
从迭代图中可以看出
- W φ W_\varphi Wφ是基于低通滤波器产生的,所以被称为approximation coefficients近似系数。
- W ψ i for i = H , V , D W_\psi^i \space \text{for }i=H,V,D Wψi for i=H,V,D分别是水平、垂直、对角的detail coefficients。
由于FWT是向创建低分辨率分量, f ( x , y ) f(x,y) f(x,y)表示高分辨率分量,所以在第一层叠戴循环时,我们可以令 W φ ( j + 1 , m , n ) = f ( x , y ) W_\varphi(j+1,m,n)=f(x,y) Wφ(j+1,m,n)=f(x,y)
2.1 FWTs Using the Wavelet Toobox
在MATLAB中,我们若想生成FWT的filter,可以使用函数wfilters:
[Lo_D, Hi_D, Lo_R, Hi_R] = wfilters(wname)
其中,Lo_D、Hi_D、Lo_R、Hi_R都是行向量,分别为低通分解、高通分解、低通重建、高通重建滤波器;wname指的是小波滤波系数,具体如下表所示:
Wavelet | wfamily | wname |
---|---|---|
Haar | ‘haar’ | ‘haar’ |
Daubechies | ‘db’ | ‘db2’、‘db3’、…‘db5’ |
Coiflets | ‘coif’ | ‘coif1’、‘coif2’、…‘coif5’ |
Symlets | ‘sym’ | ‘sym1’、‘sym2’、…‘sym5’ |
Discrete Meyer | ‘dmey’ | ‘dmey’ |
Biorthogonal | ‘bior’ | ‘bior1.1’、‘bior1.3’、‘bior1.5’、‘bior2.2’… |
Reverse Biorthogonal | ‘rbio’ | ‘rbior1.1’、‘rbior1.3’、‘rbior1.5’、‘rbior2.2’… |
想查询更多小波族相关信息,可以在MATLAB中输入命令:
waveinfo(wfamily)
Frequently coupled filter pairs频率耦合滤波器对用如下语句代替:
[F1, F2] = wfilters(wname, type)
其中,type设置为’d’、‘r’、‘l’、'h’分别得到分解、重建、低通、高通滤波器。分解或低通滤波器将返回到F1,结果放在F2中。
如果想获得归一化正交变换的尺度和小波函数的数字近似值,可以使用函数wavefun:
[phi, psi, xval] = wavefun(wname, iter)
其中,phi和psi为返回的近似向量,phi为尺度scale,psi为小波wavelet,xval为评估向量。iter比啊是迭代次数。对于双正交变换,命令则变为:
[phi1, psi1, phi2, psi2, xval] = wavefun(wname, iter)
其中,phi1和psi1是decomposition function分解函数,phi2和psi2是reconstruction function重建函数。
Example Haar filters, scaling, and wavelet functions
>> [Lo_D, Hi_D, Lo_R, Hi_R] = wfilters('haar')
Lo_D =
0.7071 0.7071
Hi_D =
-0.7071 0.7071
Lo_R =
0.7071 0.7071
Hi_R =
0.7071 -0.7071
>> [phi, psi, xval] = wavefun('haar', 10);
>> xaxis = zeros(size(xval));
>> subplot(121);plot(xval, phi, 'k', xval, xaxis, '--k');
>> axis([0 1 -1.5 1.5]); axis square
>> title('Haar Scaling Function');
>> subplot(122);plot(xval, psi, 'k', xval, xaxis, '--k');
>> axis([0 1 -1.5 1.5]); axis square
>> title('Haar Wavelet Function');
![image-20211201205126225](/Users/wanghe/Library/Application%20Support/typora-user-images/image-20211201205126225.png)
其中,Haar的尺度函数和小波函数都是不连续且紧密支撑的,支撑是1,支撑区间的优先区间之外是0。
有了wavefilter,计算wavelet transform小波变换一般使用函数wavedec2:
[C, S] = wavedec2(X, N, Lo_D, Hi_D)
[C, S] = wavedec2(X, N, wname)
其中,X是输入图像,N是被计算的尺度数,输出的C和S分别为行向量C和记录C中系数排列的矩阵S。
Example A simple FWT using Haar filters:
>> f = magic(4)
f =
16 2 3 13
5 11 10 8
9 7 6 12
4 14 15 1
>> [c1, s1] = wavedec2(f, 1, 'haar')
c1 =
Columns 1 through 9
17.0000 17.0000 17.0000 17.0000 1.0000 -1.0000 -1.0000 1.0000 4.0000
Columns 10 through 16
-4.0000 -4.0000 4.0000 10.0000 6.0000 -6.0000 -10.0000
s1 =
2 2
2 2
4 4
其中,f是一个4x4的图像,被变换为1x16的小波分解向量c1和3x2的记录矩阵s1。
记录矩阵s1记录了矩阵的大小:每个向量的第一个元素是参考细节或近似矩阵的行数,第二个元素是列数。
- s1(end, :)记录了原始图像f的大小;
- s1(2, :)记录了3个细节系数矩阵的大小;
- s1(1, :)记录了计算近似矩阵的大小;
当我们的计算尺度为两个尺度时:
>> [c2, s2] = wavedec2(f, 2, 'haar') % 分解2次
c2 =
Columns 1 through 9
34.0000 0 0 0.0000 1.0000 -1.0000 -1.0000 1.0000 4.0000
Columns 10 through 16
-4.0000 -4.0000 4.0000 10.0000 6.0000 -6.0000 -10.0000
s2 =
1 1
1 1
2 2
4 4
这里,s2记录矩阵记录的数据为:
- s1(end, :)记录了原始图像f的大小;
- s1(3, :)记录了尺度/分解级别为2时,3个细节系数矩阵的大小;
- s1(2, :)记录了尺度/分解级别为1(第一次)时,3个细节系数矩阵的大小;
- s1(1, :)记录了最后近似矩阵的大小;
2.2 FWTs without the Wavelet Toolbox
function wavefilter:创建小波滤波器
function [varargout] = wavefilter(wname, type)
%WAVEFILTER Create wavelet decomposition and reconstruction filters .
% [VARARGOUT] = WAVEFILTER (WNAME, TYPE) returns the decomposition
% and/or reconstruction filters used in the computation of the
% forward and inverse FWT(fast wavelet trans form).
%
% EXAMPLES :
% [ld, hd, lr, hr] = wavefilter('haar') Get the low and highpass
% decomposition (ld, hd)
% and reconstruction
% (lr, hr) filters for
% wavelet 'haar'.
% [ld, hd] = wavefilter('haar', 'd') Get decomposition filters
% ld and hd.
% [lr, hr] = wavefilter('haar', 'r') Get recons truction
% filters lr and hr.
%
% INPUTS:
% WNAME Wavelet Name
% -----------------------------------------------------------------
% 'haar' or 'db1' Haar
% 'db4' 4th order Daubechies
% 'sym4' 4th order Symlets
% 'bior6.8' Cohen-Daubechies- Feauveau biorthogonal
% 'jpeg9.7' Antonini-Barlaud-Mathieu-Daubechies
%
% TYPE Filter Type
% -----------------------------------------------------------------
% 'd' Decomposition filters
% 'r' Reconstruction filters
% Check the input and output arguments.
error(nargchk(1, 2, nargin));
if (nargin == 1 && nargout ~= 4) || (nargin == 2 && nargout ~= 2)
error('Invalid number of output arguments.');
end
if nargin == 1 && ~ischar(wname)
error('WNAME must be a string.');
end
if nargin == 2 && ~ischar(type)
error('TYPE must be a string.');
end
% Create filters for the requested wavelet.
switch lower(wname)
case {'haar', 'db1'}
ld = [1 1]/sqrt(2);
hd = [-1 1]/sqrt(2);
lr = ld;
hr = -hd;
case 'db4'
ld = [-1.059740178499728e-002 3.288301166698295e-002 ...
3.084138183598697e-002 -1.870348117188811e-001 ...
-2.798376941698385e-002 6.308807679295904e-001 ...
7.148465705525415e-001 2.303778133088552e-001];
t = (0:7);
hd = ld;
hd(end:-1:1) = cos(pi*t).*ld;
lr = ld;
lr(end:-1:1) = ld;
hr = cos(pi*t).*ld;
case 'sym4'
ld = [-7.576571478927333e-002 -2.963552764599851e-002 ...
4.976186676320155e-001 8.037387518059161e-001 ...
2.978577956052774e-001 -9.921954357684722e-002 ...
-1.260396726203783e-002 3.222310060404270e-002];
t = (0:7);
hd = ld;
hd(end:-1:1) = cos(pi*t).*ld;
lr = ld;
lr(end:-1:1) = ld;
hr = cos(pi*t).*ld;
case 'bior6.8'
ld = [0 1.908831736481291e-003 -1.914286129088767e-003 ...
-1.699063986760234-002 1.193456527972926e-002 ...
4.973290349094079e-002 -7.726317316720414e-002 ...
-9.405920349573646-002 4.207962846098268e-001 ...
8.259229974584023e-001 4.207962846098268e-001 ...
-9.405920349573646e-002 -7.726317316720414e-002 ...
4.973290349094079e-002 1.193456527972926e-002 ...
-1.699063986760234e-002 -1.914286129088767e-003 ...
1.908831736481291e-003];
hd = [0 0 0 1.442628250562444e-002 -1.446750489679015e-002 ...
-7.872200106262882e-002 4.036797903033992e-002 ...
4.178491091502746e-001 -7.589077294536542e-001 ...
4.178491091502746e-001 4.036797903033992e-002 ...
-7.872200106262882e-002 -1.446750489679015e-002 ...
1.442628250562444e-002 0 0 0 0];
t = (0:17);
lr = cos(pi*(t+1)).*hd;
hr = cos(pi*t).*ld;
case 'jpeg9.7'
ld = [0 0.02674875741080976 -0.01686411844287495 ...
-0.07822326652898785 0.2668641184428723 ...
0.6029490182363579 0.2668641184428723 ...
-0.07822326652898785 -0.01686411844287495 ...
0.02674875741080976];
hd = [0 -0.09127176311424948 0.05754352622849957 ...
0.5912717631142470 -1.115087052456994 ...
0.5912717631142470 0.05754352622849957 ...
-0.09127176311424948 0 0] ;
t = (0:9);
lr = cos(pi*(t+1)).*hd;
hr = cos(pi*t).*ld;
otherwise
error('Unrecognizable wavelet name (WNAME).');
end
% Output the requested filters.
if (nargin == 1)
varargout(1:4) = {ld, hd, lr, hr};
else
switch lower(type(1))
case 'd'
varargout = {ld, hd};
case 'r'
varargout = {lr, hr};
otherwise
error('Unrecognizable filter TYPE');
end
end
function wavefast:实现FWT
function [c, s] = wavefast(x, n, varargin)
%WAVEFAST Computes the FWT of a '3-D extended' 2-D array.
% [C, L] = WAVEFAST(X, N, LP, HP) computes ' PAGES' 2D N-level
% FWTs of a 'ROWS x COLUMNS x PAGES' matrix X with respect to
% decomposition filters LP and HP .
%
% [C, L] = WAVEFAST(X, N, WNAME) performs the same operation but
% fetches filters LP and HP for wavelet WNAME using WAVEFILTER.
% Scale parameter N must be less than or equal to log2 of the
% maximum image dimension. Filters LP and HP must be even. To
% reduce border distortion, X is symmetrically extended. That is,
% if X = [c1 c2 c3 ... cn] (in 1D), then its symmetric extension
% would be [... c3 c2 c1 c1 c2 c3 ... cn cn cn-1 cn-2 ...].
%
% OUTPUTS:
% Vector C is a coef f icient decomposition vector:
%
% C = [ a1(n) ... ak(n) h1(n) ... hk(n) v1(n) ... vk(n)
% d1 (n) ... dk(n) h1 (n-1) ... d1(1) ... dk(1) ]
%
% where ai, hi, vi, and di for i = 0,1, ... k are columnwise
% vectors containing approximation, horizontal, vertical, and
% diagonal coefficient matrices, respectively, and k is the
% number of pages in the 3-D extended array X. C has 3n + 1
% sections where n is the number of wavelet decomposi tions .
%
% Matrix S is an [(n+2) x 2] bookkeeping matrix if k = 1;
% else it is [(n+2) x 3]:
%
% S =[ sa(n, :); sd(n, :); sd(n-1, :); ... ; sd(1, :); sx ]
%
% where sa and sd are approximation and detail size entries.
% Check the input arguments for reasonableness
error(nargchk(3, 4, nargin));
if nargin == 3
if ischar(varargin{1})
[lp, hp] = wavefilter(varargin{1}, 'd');
else
error('Missing wavelet name');
end
else
lp = varargin{1};
hp = varargin{2};
end
% Get the filter length, 'lp', input array size, 'sx', and number of pages,
% 'pages', in extended 2-D array x
fl = length(lp);
sx = size(x);
pages = size(x, 3);
if ((ndims(x) ~= 2) && (ndims(x) ~= 3)) || (min(sx) < 2) || ~isreal(x) ...
|| ~isnumeric(x)
error('X must be a real, numeric 2-D or 3-D matrix.');
end
if (ndims(lp) ~= 2) || ~isreal(lp) || ~isnumeric(lp) ...
|| (ndims(hp) ~= 2) || ~isreal(hp) || ~isnumeric(hp) ...
|| (fl ~= length(hp)) || rem(fl, 2) ~= 0
error(['LP and HP must be even and equal length real, ' ...
'numeric filter vectors.']);
end
if ~isreal(n) || ~isnumeric(n) || (n < 1) || (n > log2(max(sx)))
error('N must be a real scaler between 1 and log2(max(size(X)))');
end
% Init the starting output data structures and initial approximation
c = [];
s = sx(1:2);
app = cell(pages, 1);
for i = 1:pages
app{i} = double(x(:, :, i));
end
% For each decomposition ...
for i = 1:n
% Extend the appriximation symmetrically.
[app, keep] = symextend(app, fl, pages);
% Convolve rows with HP and downsample. Then convolve columns with HP
% and LP to get the diagonal and vertical coefficients.
rows = symconv(app, hp, 'row', fl, keep, pages);
coefs = symconv(rows, hp, 'col', fl, keep, pages);
c = addcoefs(c, coefs, pages);
s = [size(coefs{1}); s];
coefs = symconv(rows, lp, 'col', fl, keep, pages);
c = addcoefs(c, coefs, pages);
% Convolve rows with LP and downsample. Then convolve columns with HP
% and LP to get the horizontal and next approximation coeffcients.
rows = symconv(app, lp, 'row', fl, keep, pages);
coefs = symconv(rows, hp, 'col', fl, keep, pages);
c = addcoefs(c, coefs, pages);
app = symconv(rows, lp, 'col', fl, keep, pages);
end
% Append the final appriximation structures
c = addcoefs(c, app, pages);
s = [size(app{1}); s];
if ndims(x) > 2
s(:, 3) = size(x, 3);
end
%----------------------------------------------------------------------
function nc = addcoefs(c, x, pages)
% Add 'pages' array coefficients to the wavelet decomposition vector.
nc = c;
for i = pages:-1:1
nc = [x{i}(:)' nc];
end
%----------------------------------------------------------------------
function [y, keep] = symextend(x, fl, pages)
% Compute the number of coefficients to keep after convolution and
% downsampling. Then extend the 'pages' arrays of x in both dimensions.
y = cell(pages, 1);
for i = 1:pages
keep = floor((fl + size(x{i}) - 1) / 2);
y{i} = padarray(x{i}, [(fl - 1) (fl - 1)], 'symmetric', 'both');
end
%----------------------------------------------------------------------
function y = symconv(x, h, type, fl, keep, pages)
% For the 'pages' 2-D arrays in x convolve the rows or columns with h,
% downsample, and extract the center section since symmetrically extened.
y = cell(pages, 1);
for i = 1:pages
if strcmp(type, 'row')
y{i} = conv2(x{i}, h);
y{i} = y{i}(:, 1:2:end);
y{i} = y{i}(:, fl / 2 + 1:fl / 2 + keep(2));
else
y{i} = conv2(x{i}, h');
y{i} = y{i}(1:2:end, :);
y{i} = y{i}(fl / 2 + 1:fl / 2 + keep(1), :);
end
end
function fwtcompare:自定义的FWT于wavedec2性能对比
function [ratio, maxdiff] = fwtcompare(f, n, wname)
%FWTCOMPARE Compare wavedec2 and wavefast.
% [RATIO, MAXDIFF] = FWTCOMPARE(E, N, WNAME) compares the operation of
% Wavelet Toolbox function WAVEDEC2 and custom function WAVEFAST.
%
% INPUTS:
% F Image to be trans formed.
% N Number of scales to compute.
% WNAME Wavelet to use.
%
% OUTPUTS:
% RATIO Execution time ratio (custom/toolbox)
% MAXDIFF Maximum coefficient di f ference .
% Get transform and computation time for wavedec2.
w1 = @() wavedec2(f, n, wname);
reftime = timeit(w1);
% Get transform and computation time for wavefast.
w2 = @() wavefast(f, n, wname);
t2 = timeit(w2);
% Compare the results.
ratio = t2 / reftime;
maxdiff = abs(max(w1() - w2()));
对比结果如下:
>> [ratio, maxdifference] = fwtcompare(f, 5, 'db4')
>> [ratio, maxdifference] = fwtcompare(f, 5, 'db4')
ratio =
0.8640
maxdifference =
1.6467e-12
3 Working with Wavelet Decomposition Sructures
我们知道自定义和Toolbox中的小波计算函数都会返回**[c,S]形式的结构变量。其中c**是transform coefficient vector变换系数向量,S是bookkeeping matrix记录矩阵。
这里的c我们是将多尺度二维小波变换系数集成为单一的一维向量:
c
=
[
A
N
(
:
)
′
H
N
(
:
)
′
V
N
(
:
)
′
D
N
(
:
)
′
…
H
i
(
:
)
′
V
i
(
:
)
′
D
i
(
:
)
′
…
H
1
(
:
)
′
V
1
(
:
)
′
D
1
(
:
)
′
]
c=[A_N(:)' \space\space H_N(:)' \space\space V_N(:)' \space\space D_N(:)' \dots H_i(:)'\space\space V_i(:)'\space\space D_i(:)'\space\space \dots H_1(:)'\space\space V_1(:)'\space\space D_1(:)'\space\space]
c=[AN(:)′ HN(:)′ VN(:)′ DN(:)′…Hi(:)′ Vi(:)′ Di(:)′ …H1(:)′ V1(:)′ D1(:)′ ]
其中,
A
N
A_N
AN是第N个分解等级的近似系数矩阵,
H
i
H_i
Hi、
V
i
V_i
Vi、
D
i
D_i
Di分别是级别i的水平、垂直、对角的变换系数矩阵,i=1,2,…,N。
设定一共进行N级分解,也就是说,迭代循环N次,c中应包含3N+1个元素。当i=1时,计算最高尺度系数;当i=N时,计算最低尺度系数。所以,c的系数是由低尺度到高尺度排列的。
矩阵S时(N+2)x2的记录数组,形式为:
S
=
[
s
a
N
;
s
d
N
;
s
d
N
−
1
;
…
s
d
i
;
…
s
d
1
;
s
f
]
S=[sa_N;\space\space sd_N;\space\space sd_{N-1};\dots\space\space sd_i;\dots\space\space sd_1;\space\space sf]
S=[saN; sdN; sdN−1;… sdi;… sd1; sf]
s
a
N
sa_N
saN表示第N级的近似矩阵
A
N
A_N
AN,
s
d
i
sd_i
sdi表示第i级的细节矩阵大小(这里三个方向的细节矩阵都被转换为列向量进行存储),
s
f
sf
sf表示原始图像的大小。
对于RGB彩色图像,即三位数组作为输入时,我们将其作为extended 2-D array扩展的二维数组进行处理。也就是说,矩阵由矩阵构成,(实际上很简单),表达如下:
A
N
(
:
)
′
=
[
A
N
1
(
:
)
′
A
N
2
(
:
)
′
…
A
N
k
(
:
)
′
]
H
i
(
:
)
′
=
[
H
i
1
(
:
)
′
H
i
2
(
:
)
′
…
H
i
k
(
:
)
′
]
V
i
(
:
)
′
=
[
V
i
1
(
:
)
′
V
i
2
(
:
)
′
…
V
i
k
(
:
)
′
]
D
i
(
:
)
′
=
[
D
i
1
(
:
)
′
D
i
2
(
:
)
′
…
D
i
k
(
:
)
′
]
A_N(:)'=[A_N^1(:)'\space\space A_N^2(:)' \dots A_N^k(:)']\\ H_i(:)'=[H_i^1(:)'\space\space H_i^2(:)' \dots H_i^k(:)']\\ V_i(:)'=[V_i^1(:)'\space\space V_i^2(:)' \dots V_i^k(:)']\\ D_i(:)'=[D_i^1(:)'\space\space D_i^2(:)' \dots D_i^k(:)']
AN(:)′=[AN1(:)′ AN2(:)′…ANk(:)′]Hi(:)′=[Hi1(:)′ Hi2(:)′…Hik(:)′]Vi(:)′=[Vi1(:)′ Vi2(:)′…Vik(:)′]Di(:)′=[Di1(:)′ Di2(:)′…Dik(:)′]
其中,k时扩展数组的页数,i时分解的级数。这样我们得到的**[c,S]**,c还是有3N+1个元素构成,S则变成了(N+2)x3个数组。在数组中,第三列指定了c中二位数组的数目。
Example Wavelet Toolbox functions for manipulating transform decomposition vector c:
>> f = magic(8);
>> [c1, s1] = wavedec2(f, 3, 'haar');
>> size(c1)
ans =
1 64
>> s1
s1 =
1 1
1 1
2 2
4 4
8 8
>> approx = appcoef2(c1, s1, 'haar')
approx =
260.0000
>> horizdet2 = detcoef2('h', c1, s1, 2)
horizdet2 =
1.0e-13 *
0 -0.2842
0 0
>> newc1 = wthcoef2('h', c1, s1, 2);
>> newhorizdet2 = detcoef2('h', newc1, s1, 2)
newhorizdet2 =
0 0
0 0
这里,function appcoef2是根据小波分解结构[C,S],计算N级重构系数矩阵,即获得近似矩阵:
a = appcoef2(c, s, wname)
function detcoef2是获取detail matrix,语法如下:
d = detcoef2(o, c, s, n)
其中,o可以设置为’h’、’ v’、'd’分别表示水平、垂直、对角的细节矩阵。
function wthcoef2是用于对二维信号的小波系数阈值进行处理的函数:
nc = wthcoef2(type, c, s, n, t, sorh)
type针对阈值近似系数设置为’a’,细节阈值可以传参’h’、’ v’、'd’分别表示水平、垂直、对角。n时被阈值处理后的分解级别的向量,阈值以向量t中相应的阈值为基础,如果忽略t,所有满足type和n规范的系数都将被赋值为0。
3.1 Editing Wavelet Decomposition Coefficients without the Wavelet Toolbox
function wavework:
function [varargout] = wavework(opcode, type, c, s, n, x)
%WAVEWORK is used to edit wavelet decomposition structures .
% [VARARGOUT] = WAVEWORK (OPCODE, TYPE, c, s, N, X) gets the
% coefficients specified by TYPE and N for access or modi fication
% based on OPCODE.
%
% INPUTS:
% OPCODE Operation to perform
% -----------------------------------------------------------------
% 'copy' [varargout] = Y = requested (via TYPE and N)
% coefficient matrix
% 'cut' [varargout] = [NC, Y] = New decompos ition vector
% (with requested coefficient matrix zeroed) AND
% requested coefficient matrix
% 'paste' [varargout] = [NC] = new decomposition vector with
% coefficient matrix replaced by X
%
% TYPE Coefficient category
% -----------------------------------------------------------------
% 'a' Approximation coef ficients
% 'h' Horizontal details
% 'v' Vertical details
% 'd' Diagonal details
%
% [c, S] is a wavelet toolbox decomposition structure.
% N is a decomposition level (Ignored if TYPE = 'a').
% X is a 2- or 3-D coefficient matrix for pasting.
error(nargchk(4, 6, nargin));
if (ndims(c) ~= 2) || (size(c, 1) ~= 1)
error('C must be a row vector.');
end
if (ndims(s) ~= 2) || ~isreal(s) || ~isnumeric(s) || ...
((size(s, 2) ~= 2) && (size(s, 3) ~= 3))
error('S must be a real, numeric two- or three-column array.');
end
elements = prod(s, 2);
if (length(c) < elements(end)) || ...
~(elements(1) + 3 * sum(elements(2:end - 1)) >= elements(end))
error('[C S] must form a standard wavelet decomposition structre.');
end
if strcmpi(opcode(1:3), 'pas') && nargin < 6
error('Not enough input arguments');
end
if nargin < 5
n = 1;
end
nmax = size(s, 1) - 2;
aflag = (lower(type(1)) == 'a');
if ~aflag && (n > nmax)
error('N exceeds the decompositions in [C S].');
end
switch lower(type(1)) % Make pointers into C
case 'a'
nindex = 1;
start = 1;
stop = elements(1);
ntst = nmax;
case {'h', 'v' ,'d'}
switch type
case 'h'
offset = 0; % Offset to details
case 'v'
offset = 1;
case 'd'
offset = 2;
end
nindex = size(s, 1) - n; % Index to detail infow
start = elements(1) + 3 * sum(elements(2:nmax - n + 1)) + ...
offset * elements(nindex) + 1;
stop = start + elements(nindex) - 1;
ntst = n;
otherwise
error('TYPE must begin with "a", "h", "v" or "d".');
end
switch lower(opcode) % Do requested action
case {'copy', 'cut'}
y = c(start:stop);
nc = c;
y = reshape(y, s(nindex, :));
if strcmpi(opcode(1:3), 'cut')
nc(start:stop) = 0;
varargout = {nc, y};
else
varargout = {y};
end
case 'paste'
if numel(x) ~= elements(end - ntst)
error('X is not sized for the requested paste.');
else
nc = c;
nc(start:stop) = x(:);
varargout = {nc};
end
otherwise
error('Unrecognized OPCODE.');
end
function wavecut
function [nc, y] = wavecut(type, c, s, n)
%WAVECUT Zeroes coefficients in a wavelet decomposition s tructure .
% [NC, Y] = WAVECUT (TYPE, C, s, N) returns a new decomposition
% vector whose detail or approximation coefficients (based on TYPE
% and N) have been zeroed. The coefficients that were zeroed are
% re turned in Y.
%
% INPUTS:
% TYPE Coefficient category
% -----------------------------------
% 'a' Approximation coefficients
% 'h' Horizontal details
% 'v' Vertical details
% 'd' Diagonal details
%
% [c, S] is a wavelet data structure.
% N specifies a decomposition level (ignored if TYPE = 'a') .
error(nargchk(3, 4, nargin));
if nargin == 4
[nc, y] = wavework('cut', type, c, s, n);
else
[nc, y] = wavework('cut', type, c, s);
end
function wavepaste
function nc = wavepaste(type, c, s, n, x)
%WAVEPASTE Puts coefficients in a wavelet decomposition structure .
% NC = WAVEPASTE(TYPE, c, S, N, X) returns the new decomposition
% structure after pasting X into it based on TYPE and N.
%
% INPUTS:
% TYPE Coefficient category
% ----------------------------------------
% 'a' Approximation coefficients
% 'h' Horizontal details
% 'v' Vertical details
% 'd' Diagonal details
%
% [c, S] is a wavelet data structure.
% N specifies a decomposition level (Ignored if TYPE = 'a').
% X is a 2- or 3-D approximation or detail coefficient
% matrix whose dimensions are appropriate for decomposition
% level N.
error(nargchk(5, 5, nargin));
nc = wavework('paste', type, c, s, n, x);
function wavecopy
function y = wavecopy(type, c, s, n)
%WAVECOPY Fetches coefficients of a wavelet decomposition structure.
% Y = WAVECOPY (TYPE, c, S, N) returns a coefficient array based on
% TYPE and N.
%
% INPUTS:
% TYPE Coefficient category
% ----------------------------------------
% 'a' Approximation coefficients
% 'h' Horizontal details
% 'v' Vertical details
% 'd' Diagonal details
%
% [c, S] is a wavelet data structure .
% N specifies a decomposition level (ignored if TYPE = 'a') .
error(nargchk(3, 4, nargin));
if nargin == 4
y = wavework('copy', type, c, s, n);
else
y = wavework('copy', type, c, s);
end
3.2 Displaying Wavelet Decomposition Coefficients
function wavedisplay
function w = wavedisplay(c, s, scale, border)
%WAVEDISPLAY Display wavelet decomposition coefficients.
% W = WAVEDISPLAY (c, S, SCALE, BORDER) displays and returns a
% wavelet coef ficient image.
%
% EXAMPLES:
% wavedisplay(c, s); Display w/defaults.
% foo = wavedisplay(c, s); Display and return.
% foo = wavedisp1ay(c, s, 4); Magnify the details.
% foo = wavedisplay(c, s, -4); Magnify absolute values.
% foo = wavedisplay(c, s, 1, ' append'); Keep border values.
%
% INPUTS/ OUTPUTS:
% [c, S] is a wavelet decomposition vector and bookkeeping
% matrix.
%
% SCALE Detail coefficient scaling
% -------------------------------------------------------
% 0 or 1 Maximum range (default)
% 2,3... Magnify default by the scale factor
% -1, -2... Magnify absolute values by abs (scale)
%
% BORDER Border between wavelet decompositions
% -------------------------------------------------------
% 'absorb' Border replaces image (default)
% 'append' Border increases width of image
%
% Image W: ------- ------- --------------- ------------------
% | | | |
% | a(n) | h(n) | |
% | | | |
% ------- ------- h(n-1) |
% | | | |
% | v(n) | d(n) | | h(n-2)
% | | | |
% ------- ------- ---------------
% | | |
% | v(n-1) | d(n-1) |
% | | |
% ------- ------- --------------- ------------------
% | |
% | v(n-2) | d(n-2)
% | |
%
% Here, n denotes the decomposition step scale and a, h, v, d are
% approximation, horizontal, vertical, and diagonal detail
% coefficients, respectively.
% Check input arguments for reasonableness .
error(nargchk(2, 4, nargin));
if (ndims(c) ~= 2) || (size(c, 1) ~= 1)
error('C must be a row vector.');
end
if (ndims(s) ~= 2) || ~isreal(s) || ~isnumeric(s) || ...
((size(s, 2) ~= 2) && (size(s, 2) ~= 3))
error('S must be a real, numeric, two- or three-column array.');
end
elements = prod(s, 2);
if (length(c) < elements(end)) || ...
~(elements(1) + 3 * sum(elements(2:end - 1)) >= elements(end))
error('[C S] must be a standard wavelet decomposition structure');
end
if (nargin > 2) && (~isreal(scale) || ~isnumeric(scale))
error('SCALE must be a real, numeric scalar');
end
if (nargin > 3) && (~ischar(border))
error('BORDER must be character string.');
end
if nargin == 2
scale = 1; % Default scale
end
if nargin < 4
border = 'absorb';
end
% Scale coefficients and determine pad fill
absflag = scale < 0;
scale = abs(scale);
if scale == 0
scale = 1;
end
[cd, w] = wavecut('a', c, s);
w = mat2gray(w);
cdx = max(abs(cd(:))) / scale;
if absflag
cd = mat2gray(abs(cd), [0, cdx]);
fill = 0;
else
cd = mat2gray(cd, [-cdx, cdx]);
fill = 0.5;
end
% Build gray image one decomposition at a time
for i = size(s, 1) - 2:-1:1
ws = size(w);
h = wavecopy('h', cd, s, i);
pad = ws - size(h);
frontporch = round(pad / 2);
h = padarray(h, frontporch, fill, 'pre');
h = padarray(h, pad - frontporch, fill, 'post');
v = wavecopy('v', cd, s, i);
pad = ws - size(v);
frontporch = round(pad / 2);
v = padarray(v, frontporch, fill, 'pre');
v = padarray(v, pad - frontporch, fill, 'post');
d = wavecopy('d', cd, s, i);
pad = ws - size(d);
frontporch = round(pad / 2);
d = padarray(d, frontporch, fill, 'pre');
d = padarray(d, pad - frontporch, fill, 'post');
% Add 1 pixel white border and concatenate coefficients.
switch lower(border)
case 'append'
w = padarray(w, [1 1], 1, 'post');
h = padarray(h, [1 0], 1, 'post');
v = padarray(v, [0 1], 1, 'post');
case 'absorb'
w(:, end, :) = 1;
w(end, :, :) = 1;
h(end, :, :) = 1;
v(:, end, :) = 1;
otherwise
error('Unrecongnized BORDER parameter.');
end
w = [w h; v d];
end
% Display result. If the reconstruction is an extended 2-D array
% with 2 or more pages, display as a time sequence.
if nargout == 0
if size(s, 2) == 2
imshow(w);
else
implay(w);
end
end
Example Transform coefficient display using wavedisplay:
>> f = imread('Fig0704.tif');
>> [c, s] = wavefast(f, 2, 'db4');
>> wavedisplay(c, s);
>> subplot(131);wavedisplay(c, s);
>> subplot(132);wavedisplay(c, s, 8);
>> subplot(133);wavedisplay(c, s, -8);
![image-20211202231559366](/Users/wanghe/Library/Application%20Support/typora-user-images/image-20211202231559366.png)
在不加比例进行缩放时,细节系数在第一张图中很难看出;
在第二张图中,图像被放大了8倍得以强调,这里级别为1的系数子图像边界的中等灰度填充,插入是为了调节变换系数子图的维数变换;
第三张图显示了细节绝对值的结果。
4 The Inverse Fast Wavelet Transform
类似于FWT正向变换,只不过顺序是颠倒了一下,而且,FWT正向变换中用到的downsampling下采样要转换为相应的upsampling上采样。
![image-20211202232224473](/Users/wanghe/Library/Application%20Support/typora-user-images/image-20211202232224473.png)
用公式表达如下:
I
D
=
W
ψ
D
(
j
,
m
,
n
)
↑
m
=
k
/
2
,
k
≥
0
⊙
h
ψ
(
m
)
I
V
=
W
ψ
V
(
j
,
m
,
n
)
↑
m
=
k
/
2
,
k
≥
0
⊙
h
φ
(
m
)
I
H
=
W
ψ
H
(
j
,
m
,
n
)
↑
m
=
k
/
2
,
k
≥
0
⊙
h
ψ
(
m
)
I
A
=
W
φ
(
j
,
m
,
n
)
↑
m
=
k
/
2
,
k
≥
0
⊙
h
φ
(
m
)
W
φ
(
j
+
1
,
m
,
n
)
=
[
[
I
D
+
I
V
]
↑
n
=
k
/
2
,
k
≥
0
]
⊙
h
ψ
(
n
)
+
[
[
I
H
+
I
A
]
↑
n
=
k
/
2
,
k
≥
0
]
⊙
h
φ
(
n
)
I_D=W_\psi^D(j,m,n)\uparrow^{m=k/2,k\ge0}\odot h_\psi(m)\\ I_V=W_\psi^V(j,m,n)\uparrow^{m=k/2,k\ge0}\odot h_\varphi(m)\\ I_H=W_\psi^H(j,m,n)\uparrow^{m=k/2,k\ge0}\odot h_\psi(m)\\ I_A=W_\varphi(j,m,n)\uparrow^{m=k/2,k\ge0}\odot h_\varphi(m)\\ W_\varphi(j+1,m,n)=[[I_D+I_V]\uparrow^{n=k/2,k\ge0}]\odot h_\psi(n)+[[I_H+I_A]\uparrow^{n=k/2,k\ge0}]\odot h_\varphi(n)
ID=WψD(j,m,n)↑m=k/2,k≥0⊙hψ(m)IV=WψV(j,m,n)↑m=k/2,k≥0⊙hφ(m)IH=WψH(j,m,n)↑m=k/2,k≥0⊙hψ(m)IA=Wφ(j,m,n)↑m=k/2,k≥0⊙hφ(m)Wφ(j+1,m,n)=[[ID+IV]↑n=k/2,k≥0]⊙hψ(n)+[[IH+IA]↑n=k/2,k≥0]⊙hφ(n)
公式也可简化为:
D
→
h
p
、
h
p
V
→
l
p
、
h
p
H
→
h
p
、
l
p
A
→
l
p
、
l
p
D\rightarrow hp、hp\\ V\rightarrow lp、hp\\ H\rightarrow hp、lp\\ A\rightarrow lp、lp
D→hp、hpV→lp、hpH→hp、lpA→lp、lp
一般计算IFWT时,需要使用function waverec2:
g = waverec2(C, S, wname)
其中g是得到的重构后的二维图像。要求的重构滤波器可以使用如下语句:
g = waverec2(C, S, LO_R, HI_R)
我们也可以根据IFWT的迭代图,自定义一个逆快速小波变换的函数function waveback
function [varargout] = waveback(c, s, varargin)
%WAVEBACK Computes inverse FWTs for multi-level decomposition [c, S] .
% [VARARGOUT] = WAVEBACK(C, s, VARARGIN) performs a 2D N-level
% partial or complete wavelet reconstruction of decomposition
% structure [C, S].
%
% SYNTAX:
% Y = WAVEBACK(C, S, ' WNAME'); Output inverse FWT matrix Y
% Y = WAVEBACK(C, S, LR, HR); using lowpass and highpass
% reconstruction filters (LR and
% HR) or filters obtained by
% calling WAVEFILTER with 'WNAME'.
%
% [NC, NS] = WAVEBACK(C, S, 'WNAME', N); Output new wavelet
% [NC, NS] = WAVEBACK(C, s, LR, HR, N); decomposition structure
% [NC, NS] after N step
% reconstruction.
% Check the input and output arguments for reasonableness.
error(nargchk(3, 5, nargin));
error(nargchk(1, 2, nargout));
if (ndims(c) ~= 2) || (size(c, 1) ~= 1)
error('C must be a row vector.');
end
if (ndims(s) ~= 2) || ~isreal(s) || ~isnumeric(s) || ...
((size(s, 2) ~= 2) && (size(s, 2) ~= 3))
error('S must be a real, numeric two- or three column array.');
end
elements = prod(s, 2);
if (length(c) < elements(end)) || ...
~(elements(1) + 3 * sum(elements(2:end- 1 )) >= elements(end))
error('[C, S] must be a standard wavelet decomposition structure.');
end
% Maximum levels in [C, S].
nmax = size(s, 1) - 2;
% Get third input parameter and init check flags.
wname = varargin{1};
filterchk = 0;
nchk = 0;
switch nargin
case 3
if ischar(wname)
[lp, hp] = wavefilter(wname, 'r');
n = nmax;
else
error('Undefined filter.');
end
if nargout ~= 1
error('Wrong number of output arguments.');
end
case 4
if ischar(wname)
[lp, hp] = wavefilter(wname, 'r');
n = varargin{2};
nchk = 1;
else
lp = varargin{1};
hp = varargin{2};
filterchk = 1;
n = nmax;
if nargout ~= 1
error('Wrong number of output arguments.');
end
end
case 5
lp = varargin{1};
hp = varargin{2};
filterchk = 1;
n = varargin{3};
nchk = 1;
otherwise
error('Improper number of input arguments');
end
fl = length(lp);
if filterchk
if (ndims(lp) ~= 2) || ~isreal(lp) || ~isnumeric(lp) || ...
(ndims(hp) ~= 2) || ~isreal(hp) || ~isnumeric(hp) || ...
(fl ~= length(hp)) || (rem(fl, 2) ~= 0)
error('LP and HP must be even and equal length real numeric filter vectors');
end
end
if nchk && (~isnumeric(n) || ~isreal(n))
error('N must be a real numeric');
end
if (n > nmax) || (n < 1)
error('Invalid number(N) of reconstructions requested.');
end
if (n ~= nmax) && (nargout ~= 2)
error('Not enough output arguments.');
end
% Init decomposition
nc = c;
ns = s;
nnmax = nmax;
for i = 1:n
% Compute a new approximation
a = symconvup(wavecopy('a', nc, ns), lp, lp, fl, ns(3, :)) + ...
symconvup(wavecopy('h', nc, ns, nnmax), ...
hp, lp, fl, ns(3, :)) + ...
symconvup(wavecopy('v', nc, ns, nnmax), ...
lp, hp, fl, ns(3, :)) + ...
symconvup(wavecopy('d', nc, ns, nnmax), ...
hp, hp, fl, ns(3, :));
% Update decomposition
nc = nc(4 * prod(ns(1, :)) + 1:end);
nc = [a(:)' nc];
ns = ns(3:end, :);
ns = [ns(1, :); ns];
nnmax = size(ns, 1) - 2;
end
% For complete reconstructions, reformat output as 2-D.
if nargout == 1
a = nc;
nc = repmat(0, ns(1, :));
nc(:) = a;
end
varargout{1} = nc;
if nargout == 2
varargout{2} = ns;
end
% --------------------------------------------------------------------
function w = symconvup(x, f1, f2, fln, keep)
% Upsample rows and convolve columns with f1; upsample columns and convolve
% rows with f2; then extract center assuming symmetrical extension
% Process each "page" (i.e., 3rd index) of an extended 2-D array
% separately; if 'x' is 2-D, size(x, 3) = 1.
% Preallocate w.
zi = fln - 1:fln + keep(1) - 2;
zj = fln - 1:fln + keep(2) - 2;
w = zeros(numel(zi), numel(zj), size(x, 3));
for i = 1:size(x, 3)
y = zeros([2 1] .* size(x(:, :, i)));
y(1:2:end, :) = x(:,:,i);
y = conv2(y, f1');
z = zeros([1 2] .* size(y));
z(:, 1:2:end) = y;
z = conv2(z, f2);
z = z(zi, zj);
w(:, :, i) = z;
end
在自定义一个与waverec2做比较的函数function ifwtcompare
function [ratio, maxdiff] = ifwtcompare(f, n, wname)
%IFWTCOMPARE Compare waverec2 and waveback.
% [RATIO, MAXDIFF] = IFWTCOMPARE(F, N, WNAME) compares the
% operation of Wavelet Too1box function WAVEREC2 and cus tom
% function WAVEBACK.
%
% INPUTS:
% F Image to transform and inverse trans form.
% N Number of scales to compute.
% WNAME Wavelet to use.
%
% OUTPUTS:
% RATIO Execution time ratio (custom/ toolbox).
% MAXDIFF Maximum generated image difference.
% Compute the transform and get output and computation time for waverec2.
[c1, s1] = wavedec2(f, n, wname);
w1 = @() waverec2(c1, s1, wname);
reftime = timeit(w1);
% Compute the transform and get output and computation time for waveback.
[c2, s2] = wavefast(f, n, wname);
w2 = @() waveback(c2, s2, wname);
t2 = timeit(w2);
% Compare the results.
ratio = t2 / reftime;
diff = double(w1()) - w2();
maxdiff = abs(max(diff(:)));
比较结果如下:
>> [ratio, maxdifference] = ifwtcompare(f, 5, 'db4')
ratio =
1.1269
maxdifference =
3.4106e-13
可以说差距比较小的。
5 Wavelet in Image Processing
基本过程如下:
- 计算一幅image的2-D Wavelet Transform二维小波变换
- 修改transform coefficients变换系数
- 计算Inverse transform反变换
Example Wavelet directionality and edge detection:
>> f = imread('Fig0707(a).tif');
>> subplot(221);imshow(f);
>> [c, s] = wavefast(f, 1, 'sym4');
>> subplot(222);wavedisplay(c, s, -6);
>> [nc, y] = wavecut('a', c, s);
>> subplot(223);wavedisplay(nc, s, -6);
>> edges = abs(waveback(nc, s, 'sym4'));
>> subplot(224);imshow(mat2gray(edges));
![image-20211203160116273](/Users/wanghe/Library/Application%20Support/typora-user-images/image-20211203160116273.png)
Example Wavelet-based image smoothing or blurring:
与Fourier Transform对应的Wavelet Transform是有效的图像平滑和图像模糊的处理方法。
>> subplot(321);imshow(f);
>> [c, s] = wavefast(f, 4, 'sym4');
>> subplot(322);wavedisplay(c, s, 20);
>> [c, g8] = wavezero(c, s, 1, 'sym4');
>> [c, s] = wavefast(f, 4, 'sym4');
>> [c, g8] = wavezero(c, s, 1, 'sym4');
>> subplot(323);imshow(g8);
>> [c, g8] = wavezero(c, s, 2, 'sym4');
>> subplot(324);imshow(g8);
>> [c, g8] = wavezero(c, s, 3, 'sym4');
>> subplot(325);imshow(g8);
>> [c, g8] = wavezero(c, s, 4, 'sym4');
>> subplot(326);imshow(g8);
![image-20211204151550173](/Users/wanghe/Library/Application%20Support/typora-user-images/image-20211204151550173.png)
Example Progressive reconstruction:
>> [c, s] = wavefast(f, 4, 'jpeg9.7');
>> wavedisplay(c, s, 8);
>> f = wavecopy('a', c, s);
>> figure; imshow(mat2gray(f));
>> [c, s] = waveback(c, s, 'jpeg9.7', 1); % Approximation 1
>> f = wavecopy('a', c, s);
>> figure; imshow(mat2gray(f));
>> [c, s] = waveback(c, s, 'jpeg9.7', 1); % Approximation 2
>> f = wavecopy('a', c, s);
>> figure; imshow(mat2gray(f));
>> [c, s] = waveback(c, s, 'jpeg9.7', 1); % Approximation 3
>> f = wavecopy('a', c, s);
>> figure; imshow(mat2gray(f));
>> [c, s] = waveback(c, s, 'jpeg9.7', 1); % Approximation 4
>> f = wavecopy('a', c, s);
>> figure; imshow(mat2gray(f));