【图像去噪】基于小波变换、contourlet变换、contourlet-小波变换+PCA算法实现SAR图像

2021-10-12 19:56 作者:Matlab工程师 0人读过 | 我要投稿

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子带系数 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); %提取细节部分系数 dcoef = dcoef .* (abs(dcoef) > thr); %硬阈值 C( prod(S(1, :)) + 1 : end) = dcoef; wim = waverec2(C, S, 'db1'); % 重构图像 figure,imshow(wim);title(sprintf('小波去噪(PSNR = %.2f dB)', ... PSNR(wim,im)) );axis on; %************小波去噪-完****************************************************** %***************contourlet变换去噪*************************************** %参数设置 pfilt='9-7'; % LP 分解滤波器 dfilt='pkva'; % DFB 分解滤波器 nlevs = [0,3,3,4,4,5]; % nlevs: DFB分解滤波器级数向量 % 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); %阈值判断 % 重构 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的低频系数 ch1 = detcoef2('h',C_1,S_1,1); %提取尺度1的水平方向高频系数 CV1 = detcoef2('v',C_1,S_1,1); %提取尺度1的垂直方向高频系数 cd1 = detcoef2('d',C_1,S_1,1); %提取尺度1的斜线方向高频系数 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); %阈值判断 %重构 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'); % 重构图像 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的低频系数 ch1 = detcoef2('h',C_1,S_1,1); %提取尺度1的水平方向高频系数 CV1 = detcoef2('v',C_1,S_1,1); %提取尺度1的垂直方向高频系数 cd1 = detcoef2('d',C_1,S_1,1); %提取尺度1的斜线方向高频系数 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处理 %高频部分设置阈值去噪 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); %重构 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'); % 重构图像 figure,imshow(wcim_PCA);title(sprintf('小波-contourlet+PCA去噪(PSNR = %.2f dB)', ... PSNR(wcim_PCA,im)) );axis on;