【圖像去噪】基于小波變換、contourlet變換、contourlet-小波變換+PCA算法實現(xiàn)SAR圖像
1 算法介紹
模型介紹見。
2 部分代碼
clc
clear all
close all
%讀取原始圖像
im=imread('BJ256_N_5.bmp');
im=double(im)/256;
figure,imshow(im);title('原始圖像');
n = prod(size(im));
%加噪 sigma=0.09
sigma = 0.09;
nim = im + sigma * randn(size(im));
figure,imshow(nim);title(sprintf('噪聲圖像(PSNR = %.2f dB)',PSNR(im, nim)));
%************小波去噪******************************************************
%用Donoho通用閾值公式計算閾值 x為要進行處理的圖像
% ? thr = delta * sqrt( 2 * log(n))
%計算delta
[C_1, S_1] = wavedec2(nim, 1, 'db1'); ? ? ? ? ? ? ? ? ? ? ? ? ? ? %小波分解
d = C_1( prod( S_1(1,:) ) + 2 * prod( S_1(2,:) ) + 1 : end); ? ? ? ?%HH子帶系數(shù)
delta = median( abs(d) ) / 0.6745;
thr = delta * sqrt(2*log(n)); ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? %閾值
[C, S] = wavedec2(nim, 1, 'db1'); ? ? ? ? ? ? ?%小波分解
dcoef = C( prod(S(1, :)) + 1 : end); ? ? ? ? ? ?%提取細節(jié)部分系數(shù)
dcoef = dcoef .* (abs(dcoef) > thr); ? ? ? ? ? ?%硬閾值
C( prod(S(1, :)) + 1 : end) = dcoef;
wim = waverec2(C, S, 'db1'); ? ? ? ? ? ? ? ? ? ? ?% 重構(gòu)圖像
figure,imshow(wim);title(sprintf('小波去噪(PSNR = %.2f dB)', ...
? ? ? ? ? ? ?PSNR(wim,im)) );axis on;
%************小波去噪-完******************************************************
%***************contourlet變換去噪***************************************
%參數(shù)設(shè)置
pfilt='9-7'; ? ? ? ? ? ? ? ?% ?LP 分解濾波器
dfilt='pkva'; ? ? ? ? ? ? ? % ?DFB 分解濾波器
nlevs = [0,3,3,4,4,5]; ? ? ?% ?nlevs: DFB分解濾波器級數(shù)向量
% Contourlet變換
y = pdfbdec(nim, pfilt, dfilt, nlevs); ?
[c, s] = pdfb2vec(y);
%閾值估計
nvar = pdfb_nest(size(im,1), size(im, 2), pfilt, dfilt, nlevs);
cth = 3 * sigma * sqrt(nvar);
fs = s(end, 1);
fssize = sum(prod(s(find(s(:, 1) == fs), 3:4), 2));
cth(end-fssize+1:end) = (4/3) * cth(end-fssize+1:end);
c = c .* (abs(c) > cth); ? %閾值判斷
% 重構(gòu)
y = vec2pdfb(c, s);
cim = pdfbrec(y, pfilt, dfilt);
figure,imshow(cim);title(sprintf('contourlet去噪(PSNR = %.2f dB)', ...
? ? ? ? ? ? PSNR(cim,im)) );axis on;
%**********contourlet變換去噪-完******************
%*****小波-contourlet變換去噪****************************************
[C_1, S_1] = wavedec2(nim, 1, 'db1'); ? ? ? ? ? ? ? ? ? ? ? ? ? ? %小波分解
ca1 = appcoef2(C_1,S_1,'db1',1); ?%提取尺度1的低頻系數(shù)
ch1 = detcoef2('h',C_1,S_1,1); ? ?%提取尺度1的水平方向高頻系數(shù)
CV1 = detcoef2('v',C_1,S_1,1); ? ?%提取尺度1的垂直方向高頻系數(shù)
cd1 = detcoef2('d',C_1,S_1,1); ? ? %提取尺度1的斜線方向高頻系數(shù)
xhi_dirLH = dfbdec_l(ch1, dfilt, nlevs(end)); ?%水平方向高頻contourlet變換
xhi_dirHL = dfbdec_l(CV1, dfilt, nlevs(end)); ?%垂直方向高頻contourlet變換
xhi_dirHH = dfbdec_l(cd1, dfilt, nlevs(end)); ?%斜線方向高頻contourlet變換
y = {ca1,xhi_dirLH,xhi_dirHL,xhi_dirHH};
[c, s] = pdfb2vec(y);
%閾值估計
nvar = pdfb_nest1(size(im,1), size(im, 2), pfilt, dfilt, nlevs);
cth = 3 * sigma * sqrt(nvar);
fs = s(end, 1);
fssize = sum(prod(s(find(s(:, 1) == fs), 3:4), 2));
cth(end-fssize+1:end) = (4/3) * cth(end-fssize+1:end);
c = c .* (abs(c) > cth); ? %閾值判斷
%重構(gòu)
y = vec2pdfb(c, s);
ch1_rec = dfbrec_l(y{2}, dfilt);
CV1_rec = dfbrec_l(y{3}, dfilt);
cd1_rec = dfbrec_l(y{4}, dfilt);
len = S_1(1,1)*S_1(1,2);
C_1(1:len) = ca1(1:end);
C_1(len+1:2*len) = ch1_rec(1:end);
C_1(2*len+1:3*len) = CV1_rec(1:end);
C_1(3*len+1:4*len) = cd1_rec(1:end);
wcim = waverec2(C_1, S_1, 'db1'); ? ? ? ? ? ? ? ? ? ? ?% 重構(gòu)圖像
figure,imshow(wcim);title(sprintf('小波-contourlet去噪(PSNR = %.2f dB)', ...
? ? ? ? ? ? PSNR(wcim,im)) );axis on;
? ? ? ?
? ? ? ?
%************小波-contourlet變換+PCA閾值***********************************
[C_1, S_1] = wavedec2(nim, 1, 'db1'); ? ? ? ? ? ? ? ? ? ? ? ? ? ? %小波分解
ca1 = appcoef2(C_1,S_1,'db1',1); ?%提取尺度1的低頻系數(shù)
ch1 = detcoef2('h',C_1,S_1,1); ? ?%提取尺度1的水平方向高頻系數(shù)
CV1 = detcoef2('v',C_1,S_1,1); ? ?%提取尺度1的垂直方向高頻系數(shù)
cd1 = detcoef2('d',C_1,S_1,1); ? ? %提取尺度1的斜線方向高頻系數(shù)
xhi_dirLH = dfbdec_l(ch1, dfilt, nlevs(end)); ?%水平方向高頻contourlet變換
xhi_dirHL = dfbdec_l(CV1, dfilt, nlevs(end)); ?%垂直方向高頻contourlet變換
xhi_dirHH = dfbdec_l(cd1, dfilt, nlevs(end)); ?%斜線方向高頻contourlet變換
% %PCA處理
%高頻部分設(shè)置閾值去噪
for i = 1:2^nlevs(end)
? ?%LH分量
? ?LH = cell2mat(xhi_dirLH(i));
? ?[m,n] = size(LH);
? ?for j = 1:m
? ? ? ?temp1(j,:) = LH(j,:) - mean(LH(j,:));
? ?end
? ?RLH = temp1 * temp1'/n;
? ?[evLH,edLH] = eig(RLH);
? ?yLH = evLH'*temp1;
? ?clear temp1;
? ?yLHm = mean(mean(abs(yLH)));
? ?LH = LH.*(abs(LH) > yLHm);
? ?xhi_dirLH(i) = {LH};
? ?
? ?
? ? %HL分量
? ?HL = cell2mat(xhi_dirHL(i));
? ? ? ?[m,n] = size(HL);
? ?for j = 1:m
? ? ? ?temp1(j,:) = HL(j,:) - mean(HL(j,:));
? ?end
? ?RHL = temp1 * temp1'/n;
? ?[evHL,edHL] = eig(RHL);
? ?yHL = evHL'*temp1;
? ?clear temp1;
? ?yHLm = mean(mean(abs(yHL)));
? ?HL = HL.*(abs(HL) > yHLm);
? ?xhi_dirHL(i) = {HL};
? ?
? ?%HH分量
? ?HH = cell2mat(xhi_dirHH(i)); ? ?
? ? ? ?[m,n] = size(HH);
? ?for j = 1:m
? ? ? ?temp1(j,:) = HH(j,:) - mean(HH(j,:));
? ?end
? ?RHH = temp1 * temp1'/n;
? ?[evHH,edHH] = eig(RHH);
? ?yHH = evHH'*temp1;
? ?clear temp1;
? ?yHHm = mean(mean(abs(yHH)));
? ?HH = HH.*(abs(HH) > yHHm);
? ?xhi_dirHH(i) = {HH};
end
y = {ca1,xhi_dirLH,xhi_dirHL,xhi_dirHH};
[c, s] = pdfb2vec(y);
%重構(gòu)
y = vec2pdfb(c, s);
ch1_rec = dfbrec_l(y{2}, dfilt);
CV1_rec = dfbrec_l(y{3}, dfilt);
cd1_rec = dfbrec_l(y{4}, dfilt);
len = S_1(1,1)*S_1(1,2);
C_1(1:len) = ca1(1:end);
C_1(len+1:2*len) = ch1_rec(1:end);
C_1(2*len+1:3*len) = CV1_rec(1:end);
C_1(3*len+1:4*len) = cd1_rec(1:end);
wcim_PCA = waverec2(C_1, S_1, 'db1'); ? ? ? ? ? ? ? ? ? ? ?% 重構(gòu)圖像
figure,imshow(wcim_PCA);title(sprintf('小波-contourlet+PCA去噪(PSNR = %.2f dB)', ...
? ? ? ? ? ? PSNR(wcim_PCA,im)) );axis on;
3 仿真結(jié)果
?
4 參考文獻
[1]盧曉光, 韓萍, 吳仁彪,等. 基于二維小波變換和獨立分量分析的SAR圖像去噪方法[J]. 電子與信息學(xué)報, 2008.
[1]田福苓. 基于Contourlet變換域統(tǒng)計模型的SAR圖像去噪[D]. 西安電子科技大學(xué).
5 代碼下載
