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

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;