本文代碼的實現嚴重依賴前面的兩篇文章:
[一維信號的小波閾值去噪](http://blog.csdn.net/ebowtang/article/details/40481393#t14)
[小波變換一維Mallat算法的C++實現](http://blog.csdn.net/ebowtang/article/details/40433861#t0)
圖像在獲取或傳輸過程中會因各種噪聲的干擾使質量下降,這將對后續圖像的處理產生不利影響.所以必須對圖像進行去噪處理,而去噪所要達到的目的就是在較好去除噪聲的基礎上,良好的保持圖像的邊緣等重要細節.在圖像去噪領域得到廣泛的應用.本博文根據小波的分解與重構原理,實現了基于硬閾值和軟閾值函數的小波閾值去噪的C++版本,最終結果與matlab庫函數運算結果完全一致。
注本文的大部分文字提取于參考論文
# 一,小波閾值去噪基本理論
###小波閾值處理
小波閾值收縮法是Donoho和Johnstone提出的,其主要理論依據是,小波變換具有很強的去數據相關性,它能夠使信號的能量在小波域集中在一些大的小波系數中;而噪聲的能量卻分布于整個小波域內.因此,經小波分解后,信號的小波系數幅值要大于噪聲的系數幅值.可以認為,幅值比較大的小波系數一般以信號為主,而幅值比較小的系數在很大程度上是噪聲.于是,采用閾值的辦法可以把信號系數保留,而使大部分噪聲系數減小至零.小波閾值收縮法去噪的具體處理過程為:將含噪信號在各尺度上進行小波分解,設定一個閾值,幅值低于該閾值的小波系數置為0,高于該閾值的小波系數或者完全保留,或者做相應的“收縮(shrinkage)”處理.最后將處理后獲得的小波系數用逆小波變換進行重構,得到去噪后的圖像.
###閾值函數的選取
? 閾值去噪中,閾值函數體現了對超過和低于閾值的小波系數不同處理策略,是閾值去噪中關鍵的一步。設w表示小波系數,T為給定閾值,sign(*)為符號函數,常見的閾值函數有:
硬閾值函數: ? ? (小波系數的絕對值低于閾值的置零,高于的保留不變)
? ? ?
軟閾值函數: ??
? ? ??
值得注意的是:
1) 硬閾值函數在閾值點是不連續的,在下圖中已經用黑線標出。不連續會帶來振鈴,偽吉布斯效應等。
2) 軟閾值函數,原系數和分解得到的小波系數總存在著恒定的偏差,這將影響重構的精度使得重構圖像的邊緣模糊等現象.
同時這兩種函數不能表達出分解后系數的能量分布。見下圖:

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?圖1
于是不少文章出現了折衷方案,一種新的閾值函數(非本文重點):

參考代碼:
~~~
%Generate signal and set threshold.
y = linspace(-1,1,100);
subplot(311);
plot(y);title('原始線');
thr = 0.5;
% Perform hard thresholding.
ythard = wthresh(y,'h',thr);
subplot(312);
plot(ythard);title('硬閾值線');
% Perform soft thresholding.
ytsoft = wthresh(y,'s',thr);
subplot(313);
plot(ytsoft);title('軟閾值線');
~~~
### 閾值的確定
選取的閾值最好剛好大于噪聲的最大水平,可以證明的是噪聲的最大限度以非常高的概率低于
(此閾值是Donoho提出的),其中根號右邊的這個參數(叫做sigma)就是估計出來的噪聲標準偏差(第一級分解出的小波細節系數,即整個HH系數的絕對值的中值),本文將用此閾值去處理各尺度上的細節系數。
### 閾值策略
曾經做的ppt挪用過來的。

###二維信號分解與重構
小波的分解

小波的重構

# 二,Matlab庫函數實現
### 1,核心庫函數說明
1)wavedec2
圖像的多級小波分解,將返回分解出來的小波系數以及小波系數的各級長度
2)waverec2
多級小波系數的重構,重構出原信號
3)wthresh函數
對系數進行指定類型(全局閾值或者分層閾值)的閾值去噪處理
更具體的函數說明可以在,matlab里鍵入“doc 函數名”將得到很詳細的說明,當然也可以百度
### 2,軟,硬閾值處理效果:




局部放大圖像:
四幅圖象均放大兩倍,便于查看區別




### 3,完整的代碼實現
說明:代碼實現對圖像一層分解后的細節系數進行全局閾值處理(即LHL,LH,HH均用同一閾值處理),并且閾值是自定義的。
~~~
clc;
clear all;
close all;
%%%%%%%%%%%%%%%%%%%%%%%%通過matlab的函數來實現閾值去噪%%%%%%%%%%%%%%%%%%%%%%%%%% %
X=imread('sample1.bmp');%有噪聲的圖像
X=double(X);
figure(1);
imshow(X,[]);title('原始圖像'); %注意,原始圖像帶有噪聲
% 獲取輸入參數
w = 'db3';%小波類型
n = 1;%分解層數
thr = 23.91;%自定義閾值
sorh1 = 'h';%硬閾值
sorh2 = 's';%軟閾值
% 對圖像進行小波分解
[c,l] = wavedec2(X,n,w);
% 對小波系數全局閾值處理
cxch = c;% 保留近似系數
cxcs = c;% 保留近似系數
justdet = prod(l(1,:))+1:length(c);%截取細節系數(不處理近似系數)
% 閾值處理細節系數
cxch(justdet) = wthresh(cxch(justdet),sorh1,thr);%硬閾值去噪
cxcs(justdet) = wthresh(cxcs(justdet),sorh2,thr);%軟閾值去噪
%小波重建
xch = waverec2(cxch,l,w);
xcs = waverec2(cxcs,l,w);
figure(2);
imshow(uint8(xch));title('硬閾值去噪圖像');
figure(3);
imshow(uint8(xcs));title('軟閾值去噪圖像');
~~~
# 三,C加加實現
說明:如同一維的閾值去噪一樣,在執行自己編寫的wavedec2函數時必須先初始化,初始化的目的是為了獲取信號的長度,選擇的是什么小波,以及分解的等級等信息,然后計算出未來的各種信息,比如每個等級的系數的size,其中共有變量m_msgCL2D記錄了這些信息。二維小波分解的初始化函數如下:
~~~
//初始化二維圖像的分解信息,保存未來需要的信息
bool CWavelet::InitDecInfo2D(
const int height,//預分解的圖像的高度
const int width,//預分解的圖像的寬度
const int Scale,//分解尺度
const int dbn//db濾波器編號,有默認值
)
{
if (dbn != 3)
SetFilter(dbn);
if (height < m_dbFilter.filterLen - 1 || width < m_dbFilter.filterLen - 1)
{
cerr << "錯誤信息:濾波器長度大于信號的高度或者寬度!" << endl;
return false;
}
int srcHeight = height;
int srcWidth = width;
m_msgCL2D.dbn = dbn;
m_msgCL2D.Scale = Scale;
m_msgCL2D.msgHeight.resize(Scale + 2);
m_msgCL2D.msgWidth.resize(Scale + 2);
//源圖像的尺寸
m_msgCL2D.msgHeight[0] = height;
m_msgCL2D.msgWidth[0] = width;
//每一尺度上的尺寸
for (int i = 1; i <= Scale; i++)//注意:每個尺度的四個分量的長寬是一樣的
{
int exHeight = (srcHeight + m_dbFilter.filterLen - 1) / 2;//對稱拓延后系數的長度的一半
srcHeight = exHeight;
m_msgCL2D.msgHeight[i] = srcHeight;
int exWidth = (srcWidth + m_dbFilter.filterLen - 1) / 2;//對稱拓延后系數的長度一半
srcWidth = exWidth;
m_msgCL2D.msgWidth[i] = srcWidth;
}
m_msgCL2D.msgHeight[Scale + 1] = srcHeight;
m_msgCL2D.msgWidth[Scale + 1] = srcWidth;
//計算總的數據個數
int tmpAllSize = 0;
int curPartSize = 0;
int prePartSize = 0;
for (int i = 1; i <= Scale; i++)
{
curPartSize = m_msgCL2D.msgHeight[i] * m_msgCL2D.msgWidth[i];
tmpAllSize += curPartSize * 4 - prePartSize;
prePartSize = curPartSize;
}
m_msgCL2D.allSize = tmpAllSize;
m_bInitFlag2D = true;
return true;
}
~~~
### 2,核心函數的實現
### 1)二維信號的單次分解
說明:本函數建立在一維的小波分解函數基礎上(DWT)
~~~
// 二維數據的小波分解
void CWavelet::DWT2(
double *pSrcImage,//源圖像數據(存儲成一維數據,行優先存儲)
int height,//圖像的高
int width,//圖像的寬
double *pDstCeof//分解出來的圖像系數
)
{
if (!m_bInitFlag2D)
{
cerr << "錯誤信息:未初始化,無法對信號進行分解!" << endl;
return;
}
if (pSrcImage == NULL || pDstCeof == NULL)
{
cerr << "錯誤信息:dwt2數據無內存" << endl;
Sleep(3000);
exit(1);
}
int exwidth = (width + m_dbFilter.filterLen - 1) / 2 * 2;//pImagCeof的寬度
int exheight = (height + m_dbFilter.filterLen - 1) / 2 * 2;//pImagCeof的高度
double *tempImage = new double[exwidth*height];
// 對每一行進行行變換
double *tempAhang = new double[width];
double *tempExhang = new double[exwidth]; // 臨時存放每一行的處理數據
for (int i = 0; i < height; i++)
{
for (int j = 0; j < width; j++)
tempAhang[j] = pSrcImage[i*width + j];//提取每一行的數據
DWT(tempAhang, width, tempExhang);
for (int j = 0; j < exwidth; j++)
tempImage[i*exwidth + j] = tempExhang[j];
}
// 對每一列進行列變換
double *tempAlie = new double[height]; // 臨時存放每一列的轉置數據
double *tempexlie = new double[exheight]; // 臨時存放每一列的處理數據
for (int i = 0; i < exwidth; i++)
{
// 列轉置
for (int j = 0; j < height; j++)
tempAlie[j] = tempImage[j*exwidth + i];//提取每一列數據
//執行變換
DWT(tempAlie, height, tempexlie);
// 反轉置
for (int j = 0; j < exheight; j++)
pDstCeof[j*exwidth + i] = tempexlie[j];
}
AdjustData(pDstCeof, exheight, exwidth);//調整數據
delete[] tempAlie;
tempAlie = NULL;
delete[] tempexlie;
tempexlie = NULL;
delete[] tempAhang;
tempAhang = NULL;
delete[] tempExhang;
tempExhang = NULL;
delete[] tempImage;
tempImage = NULL;
}
~~~
### 2)二維信號的單次重構
說明:
~~~
//二維小波反變換
void CWavelet::IDWT2(
double *pSrcCeof, //二維源圖像系數數據
int dstHeight,//重構出來后數據的高度
int dstWidth,//重構出來后數據的寬度
double *pDstImage//重構出來的圖像
)
{
int srcHeight = (dstHeight + m_dbFilter.filterLen - 1) / 2 * 2;
int srcWidth = (dstWidth + m_dbFilter.filterLen - 1) / 2 * 2;//pSrcCeof的高度
IAdjustData(pSrcCeof, srcHeight, srcWidth);//調整成LL,HL,LH,HH
double *tempAline = new double[srcHeight]; // 臨時存放每一列的數據
double *tempdstline = new double[dstHeight]; // 臨時存放每一列的重構結果
double *pTmpImage = new double[srcWidth*dstHeight];
// 列重構
for (int i = 0; i < srcWidth; i++)//每一列
{
// 列轉置
for (int j = 0; j<srcHeight; j++)
tempAline[j] = pSrcCeof[j*srcWidth + i];//提取每一列
IDWT(tempAline, dstHeight, tempdstline);
// 反轉置
for (int j = 0; j < dstHeight; j++)
pTmpImage[j*srcWidth + i] = tempdstline[j];
}
// 對每一行進行行變換
double *tempAhang = new double[srcWidth];
double *tempdsthang = new double[dstWidth]; // 臨時存放每一行的處理數據
for (int i = 0; i < dstHeight; i++)
{
for (int j = 0; j < srcWidth; j++)
tempAhang[j] = pTmpImage[i*srcWidth + j];//提取每一行的數據
IDWT(tempAhang, dstWidth, tempdsthang);
for (int j = 0; j < dstWidth; j++)
pDstImage[i*dstWidth + j] = tempdsthang[j];
}
delete[] tempAline;
tempAline = NULL;
delete[] tempdstline;
tempdstline = NULL;
delete[] tempAhang;
tempAhang = NULL;
delete[] tempdsthang;
tempdsthang = NULL;
delete[] pTmpImage;
pTmpImage = NULL;
}
~~~
### 3)二維信號的多級分解
說明:對于每一級分解都將調用單次二維分解函數來實現,所以本函數是建立在函數IDW2基礎上
~~~
// 二維小波多級分解,需要先初始化獲取未來數據信息
bool CWavelet::WaveDec2(
double *pSrcData,//源圖像數據,存儲為一維信號
double *pDstCeof//分解后的系數,它的大小必須是m_msgCL2D.allSize
)
{
if (!m_bInitFlag2D)
{
cerr << "錯誤信息:未初始化,無法對圖像進行分解!" << endl;
return false;
}
if (pSrcData == NULL || pDstCeof == NULL)//錯誤:無內存
return false;
int height = m_msgCL2D.msgHeight[0];
int width = m_msgCL2D.msgWidth[0];
int scale = m_msgCL2D.Scale;
// 臨時變量,圖像數據
double *tempImage = new double[height*width];
int maxCoefSize =4 * m_msgCL2D.msgHeight[1] * m_msgCL2D.msgWidth[1];
double *tempDst = new double[maxCoefSize];
for (int i = 0; i < height*width; i++)
tempImage[i] = pSrcData[i];
int gap = m_msgCL2D.allSize - maxCoefSize;
for (int i = 1; i <= scale; i++)
{
DWT2(tempImage, height, width, tempDst);
// 低頻子圖像的高和寬
height = m_msgCL2D.msgHeight[i];
width = m_msgCL2D.msgWidth[i];
for (int j = 0; j < height*width; j++)
tempImage[j] = tempDst[j];//提取低頻系數(近似系數)
//
for (int j = 0, k = gap; j < 4 * height*width; j++, k++)
pDstCeof[k] = tempDst[j];//所有系數
gap -= 4 * m_msgCL2D.msgWidth[i + 1] * m_msgCL2D.msgHeight[i + 1] - height*width;
}
delete[] tempDst;
tempDst = NULL;
delete[] tempImage;
tempImage = NULL;
return true;
}
~~~
### 4)多級分解系數的重構
~~~
// 根據多級分解系數重構出二維信號,必須先初始化獲取分解信息
bool CWavelet::WaveRec2(
double *pSrcCoef,//多級分解出的源系數
double *pDstData//重構出來的信號
)
{
if (!m_bInitFlag2D)
{
cerr << "錯誤信息:未初始化,無法對信號進行分解!" << endl;
return false;
}
if (pSrcCoef == NULL || pDstData == NULL)//錯誤:無內存
return false;
int height = m_msgCL2D.msgHeight[0];
int width = m_msgCL2D.msgWidth[0];
int decLevel = m_msgCL2D.Scale;
int maxCeofSize = 4 * m_msgCL2D.msgHeight[1] * m_msgCL2D.msgWidth[1];
double *pTmpImage = new double[maxCeofSize];
int minCeofSize = 4 * m_msgCL2D.msgHeight[decLevel] * m_msgCL2D.msgWidth[decLevel];
for (int i = 0; i < minCeofSize; i++)
pTmpImage[i] = pSrcCoef[i];
int gap = minCeofSize;
for (int i = decLevel; i >= 1; i--)
{
int nextheight = m_msgCL2D.msgHeight[i - 1];//重構出來的高度
int nextwidth = m_msgCL2D.msgWidth[i - 1];//重構出來的寬度
IDWT2(pTmpImage, nextheight, nextwidth, pDstData);
if (i > 1)//i==1已經重構出來了,不再需要提取系數
{
for (int j = 0; j < nextheight*nextwidth; j++)
pTmpImage[j] = pDstData[j];
for (int j = 0; j < 3 * nextheight*nextwidth; j++)
pTmpImage[nextheight*nextwidth + j] = pSrcCoef[gap + j];
gap += 3 * nextheight*nextwidth;
}
}
delete[] pTmpImage;
pTmpImage = NULL;
return true;
}
~~~
### 3,函數正確性驗證
### 1)二維單次分解與重構測試

附帶本次測試代碼:
~~~
//小波函數的二維分解與重構
int main()
{
system("color 0A");
CWavelet cw;
double s[36] = { 1, 12, 30, 4, 5, 61, 2, 3, 41, 5, 6, 27, 3, 4, 15, 6, 72, 8, 41, 5, 6, 7, 8, 9, 5, 64, 7, 8, 9, 14, 6, 27, 8, 9, 40, 31 };
int height = 6;
int width = 6;
for (int j = 0; j < 36; j++)
{
cout << s[j] << " ";
if ((j + 1) % 6 == 0)
cout << endl;
}
cout << endl;
cout << endl;
cw.InitDecInfo2D(height, width, 1, 3);
double *dst = new double[cw.m_msgCL2D.allSize];
cw.DWT2(s, height, width, dst);
for (int j = 0; j < cw.m_msgCL2D.allSize; j++)
{
cout << dst[j] << " ";
if ((j + 1) % 10 == 0)
cout << endl;
}
cout << endl;
double *dsts = new double[36];
cw.IDWT2(dst, height, width, dsts);
for (int j = 0; j < 36; j++)
{
cout << dsts[j] << " ";
if ((j + 1) % 6 == 0)
cout << endl;
}
delete[] dst;
dst = NULL;
delete[] dsts;
dsts = NULL;
system("pasue");
return 0;
}
~~~
### 2)二維多級分解與重構測試
說明:對二維數據進行了5層分解,選取的是小波族db3

附帶本次測試代碼:
~~~
//小波函數二維多級分解與重構測試
int main()
{
system("color 0A");
CWavelet cw;
double s[48] = { 1, 12, 30, 4, 5, 61, 2, 3, 41, 5, 6, 27, 3, 4, 15, 6, 72, 8, 41, 5, 6, 7, 8, 9, 5, 64, 7, 8, 9, 14, 6, 27, 8, 9, 40, 31 ,
1, 12, 30, 4, 5, 61, 2, 3, 41, 5, 6, 27
};
int height = 6;
int width = 8;
for (int j = 0; j < 48; j++)
{
cout << s[j] << " ";
if ((j + 1) % 8 == 0)
cout << endl;
}
cout << endl;
int Scale = 5;
int dbn = 2;
cw.InitDecInfo2D(height, width, Scale, dbn);
double *dstcoef = new double[cw.m_msgCL2D.allSize];
cw.WaveDec2(s,dstcoef);
for (int i = 0; i < cw.m_msgCL2D.allSize; i++)
{
cout << dstcoef[i] << " ";
if ((i + 1) % 10 == 0)
cout << endl;
}
double *dst = new double[48];
for (int i = 0; i < 48; i++)
dst[i] = 0.0;
cw.WaveRec2(dstcoef, dst);
cout << endl; cout << endl;
for (int i = 0; i < 48; i++)
{
cout << dst[i] << " ";
if ((i + 1) % 8 == 0)
cout << endl;
}
cout << endl;
delete[] dst;
dst = NULL;
delete[] dstcoef;
dstcoef = NULL;
system("pause");
return 0;
}
~~~
### 3,閾值去噪結果:
說明:本測試只是模擬測試,對圖像的處理也是一樣的(完全一致)
### 硬閾值去噪結果

### 軟閾值去噪結果

實際的圖像處理結果為:
源噪聲圖像為:

注意以下是采用:db6,3層分解,軟閾值去噪,閾值是在前文提及的閾值基礎上縮小2.5倍得到的效果:

注意以下是采用:db6,3層分解,硬閾值去噪,閾值是在前文提及的閾值基礎上縮小2.5倍得到的效果:

附帶上述C++測試驗證主函數
~~~
//二維閾值去噪測試
int main()
{
system("color 0A");
CWavelet cw;
double s[48] = { 10, 12, 30, 4, 5, 61, 2, 3, 41, 5, 6, 27, 3, 4, 15, 6, 72, 8, 41, 5, 6, 7, 8, 9, 5, 64, 7, 8, 9, 14, 6, 27, 8, 9, 40, 31,
10, 12, 30, 4, 50, 61, 2, 3, 41, 5, 6, 27};
int height = 6;
int width = 8;
for (int j = 0; j < 48; j++)
{
cout << s[j] << " ";
if ((j + 1) % 8 == 0)
cout << endl;
}
cout << endl;
int Scale = 3;
int dbn = 3;
cw.InitDecInfo2D(height, width, Scale, dbn);
double *dstcoef = new double[48];
if (!cw.thrDenoise2D(s, dstcoef))
cerr << "Error" << endl;
for (int j = 0; j < 48; j++)
{
cout << dstcoef[j] << " ";
if ((j + 1) % 8 == 0)
cout << endl;
}
delete[] dstcoef;
dstcoef = NULL;
system("pause");
return 0;
}
~~~
附帶上述matlab驗證程序
~~~
clc;
clear all;
close all;
% %%%%%%%%%%%%%%%%%%%%%%%%通過matlab的函數來實現閾值去噪%%%%%%%%%%%%%%%%%%%%%%%%%% %
X=[ 10, 12, 30, 4, 5, 61, 2, 3;
41, 5, 6, 27, 3, 4, 15, 6;
72, 8, 41, 5, 6, 7, 8, 9;
5, 64, 7, 8, 9, 14, 6, 27;
8, 9, 40, 31,10, 12, 30, 4;
50, 61, 2, 3, 41, 5, 6, 27];
X=double(X);
% 獲取輸入參數
wname = 'db3';%小波類型
n = 3;%分解層數
sorh1 = 'h';%硬閾值
sorh2 = 's';%軟閾值
% 對圖像進行小波分解
[c,l] = wavedec2(X,n,wname);
%求取閾值
N = numel(X);
[chd1,cvd1,cdd1] = detcoef2('all',c,l,1);
cvd1=cvd1(:)';
sigma = median(abs(cvd1))/0.6745;%提取細節系數求中值并除以0.6745
thr = sigma*sqrt(2*log(N));
% 對小波系數全局閾值處理
cxch = c;% 保留近似系數
cxcs = c;% 保留近似系數
justdet = prod(l(1,:))+1:length(c);%截取細節系數(不處理近似系數)
% 閾值處理細節系數
cxch(justdet) = wthresh(cxch(justdet),sorh1,thr);%硬閾值去噪
cxcs(justdet) = wthresh(cxcs(justdet),sorh2,thr);%軟閾值去噪
%小波重建
xch = waverec2(cxch,l,wname);
xcs = waverec2(cxcs,l,wname);
~~~
注:本博文為[EbowTang](http://my.csdn.net/EbowTang)原創,后續可能繼續更新本文。如果轉載,請務必復制本條信息!
原文地址:http://blog.csdn.net/ebowtang/article/details/40481539
原作者博客:http://blog.csdn.net/ebowtang
# 參考資源:
【1】《數字圖像處理》(岡薩雷斯matlab第二版)
【2】http://ivm.sjtu.edu.cn/files/wavelet/第3章wavelet_original.pdf
【3】http://media.cs.tsinghua.edu.cn/~ahz/digitalimageprocess/chapter12/chapt12_ahz.htm#a1
【4】http://wenku.baidu.com/link?url=OYRL2n-cYkZ2J10zaMscZQ-lhR05kysQ_CaB1YM1e_aqr3DakexZRm8rtBYOHlDmxC0cNAtiCopjyog_yOIH1zliUmyz2fKfOzFTAQ1wWj3
【5】《維基百科》
【6】來自中國知網若干論文
【7】小波分析及其應用__孫延奎
【8】楊建國.小波分析及其工程應用[M].北京:機械工業出版社.2005
【9】毛艷輝.小波去噪在語音識別預處理中的應用.上海交通大學碩士學位論文.2010
【10】matlab各種函數說明,及其內部函數實現