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

最新推荐文章于 2025-12-07 21:20:15 发布

原创最新推荐文章于 2025-12-07 21:20:15 发布 · 229 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#算法 #人工智能 #机器学习 #计算机视觉

1 算法介绍

1.1 小波变换

图像的二维离散小波分解和重构过程如下图所示，分解过程可描述为：首先对图像的每一行进行 1D-DWT，获得原始图像在水平方向上的低频分量 L 和高频分量 H，然后对变换所得数据的每一列进行 1D-DWT，获得原始图像在水平和垂直方向上的低频分量 LL、水平方向上的低频和垂直方向上的高频 LH、水平方向上的高频和垂直方向上的低频 HL 以及水平和垂直方向上的的高频分量 HH。

重构过程可描述为：首先对变换结果的每一列进行以为离散小波逆变换，再对变换所得数据的每一行进行一维离散小波逆变换，即可获得重构图像。由上述过程可以看出，图像的小波分解是一个将信号按照低频和有向高频进行分离的过程，分解过程中还可以根据需要对得到的 LL 分量进行进一步的小波分解，直至达到要求。

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

对于二维图像Haar变换不再从一个方向进行滤波，而是从水平和竖直两个方向进行低通和高通滤波(水平和竖直先后不影响)，用图像表述如图所示：图中a表示原图，图b表示经过一级小波变换的结果，h1 表示水平反向的细节，v1 表示竖直方向的细节，c1表示对角线方向的细节，b表示下2采样的图像。图c中表示继续进行Haar小波变换。一级Haar小波变换实际效果如图3所示

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

1.2 contourlet变换

对信号的稀疏表示是许多信号处理及应用的基础，2004年Minh N Do、Martin Vetterli提出了一种能够较好表示二维信号的数学工具--Contourlet变换。Contourlet是用金字塔方向滤波器组(PDFB)来将图像分解成不同尺度下的方向子带的。根据PDFB的结构，PDFB是一个拉普拉斯金字塔滤波器Laplacian Pyramid (LP)和一个方向滤波器组的叠加。实验证明，Contourlet变换在图像降噪，纹理，形状的特征提取方面的性能比2-D离散小波变换有了明显的提高。

为了获得平移不变性，本章所用的Nonsubsampled Contourlet变换(NSCT)是基于Nonsubsampled金字塔(NSP)和Nonsubsampled方向滤波器(NSDFB)的一种变换。首先由NSP对输入图像进行塔形分解，分解为高通和低通两个部分，然后由NSDFB将高频子带分解为多个方向子带，低频部分继续进行如上分解。NSCT是一种新型平移不变，多尺度，多方向性的快速变换。

1. Nonsubsampled Pyramid(NSP)：

Nonsubsampled Pyramid(NSP)和Contourlet的Laplacian Pyramid(LP)多尺度分析特性不同。图像通过Nonsubsampled Pyramid(NSP)进行多尺度分解，NSP去除了上采样和下采样，减少了采样在滤波器中的失真，获得了平移不变性。NSP为具有平移不变性滤波结构的NSCT多尺度分析，可以得到与LP分解一样的多尺度分析特性。图2.4(a)处分为3个尺度。

2. Nonsubsampled方向滤波器(NSDFB)

Nonsubsampled方向滤波器(NSDFB)是一个双通道的滤波器，将分布在同方向的奇异点合成NSCT的系数。方向滤波器(DFB)是Bamberger and Smith提出的。其通过一个l层的树状结构的分解，有效的将信号分成了个子带，其频带分割成为锲形。Nonsubsampled DFB(NSDFB)为非采样，减少了采样在滤波器中的失真，获得了平移不变性。并且每个尺度下的方向子图的的大小都和原图同样大小，Contourlet变换为所有子带之和等于原图。NSCT有更多的细节得以保留，变换系数是冗余的。图2.4(b)为三个尺度下对图像频域的分割图，其中每个尺度的方向子带数目以2倍递增，以在1，2，3尺度下的方向子带数目分别为2，4，8个。

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

1.3 PCA算法

PCA的算法思路主要是：数据从原来的坐标系转换到新的坐标系，由数据本身决定。转换坐标系时，以方差最大的方向作为坐标轴方向，因为数据的最大方差给出了数据的最重要的信息。第一个新坐标轴选择的是原始数据中方差最大的方向，第二个新坐标轴选择的是与第一个新坐标轴正交且方差次大的方向。重复该过程，重复次数为原始数据的特征维数。

通过这种方式获得的新的坐标系，我们发现，大部分方差都包含在前面几个坐标轴中，后面的坐标轴所含的方差几乎为0。于是，我们可以忽略余下的坐标轴，只保留前面的几个含有绝大部分方差的坐标轴。事实上，这样也就相当于只保留包含绝大部分方差的维度特征，而忽略包含方差几乎为0的特征维度，也就实现了对数据特征的降维处理。

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;