多尺度双边滤波及基于小波变换的非线性扩散

主题：对车载手势进行基于小波变换和双边滤波的图像去噪

---

一、原理

(1)基于小波变换的非线性扩散

小波分解→高频子带非线性扩散

小波分解将图像分解为fa (低频部分)fh fv fd(包含细节信息与噪声信息的高频部分)

fh(水平高频子带):进行从左至右与从右至左的非线性扩散列操作；

fv(垂直高频子带): 进行从上至下与从下至上的非线性扩散行操作

fd(对角高频子带): 进行从左至右与从右至左的非线性扩散列操作以及进行从上至下与从下至上的非线性扩散行操作；

非线性扩散方程

(2)多尺度双边滤波

①利用高斯核函数定义空间临近度与灰度值相似度

②领域空间内与像素(x,y)相似的像素集

 

σ₀ 为灰度阈值，N为像素（x,y）的领域空间

③像素集中像素的个数与领域空间之比t

t=|Q|/|N|

设立一个比例阈值t₀

若ｔ≤t₀，则在计算过程中时,只有N中与其相似的像素参加；

若ｔ＞t₀，则在计算过程中，Ｎ中所有像素都参与运算

④双边滤波计算公式

二 具体实现

1.读图并灰度化处理

clc;
close all;
%输入原图
srcI=imread('1.jpg');
figure;
imshow(srcI);
title('原始图片');
%判断通道数并决定是否进行灰度化
srcIsize=size(srcI);
if numel(srcIsize)>2
    srcI=rgb2gray(srcI);
end
srcI=im2double(srcI);
figure;
imshow(srcI);
title('灰度化后的图片');

2.哈尔平稳小波一尺度的分解

%小波分解
[cA,cH,cV,cD]=swt2(srcI,1,'haar');
figure;
subplot(2,2,1);
imshow(cA);
title('低频图像');
subplot(2,2,2);
imshow(cH);
title('水平高频子带图像');
subplot(2,2,3);
imshow(cV);
title('垂直高频子带图像');
subplot(2,2,4);
imshow(cD);
title('对角高频子带图像');

3.不同方向的非线性扩散

%非线性扩散
k=15;
lamda=0.15;
t=20;
%对水平高频子带进行列操作
dstCH=columchange(cH,k,lamda,t);
%对垂直高频子带进行行操作
dstCV=rowchange(cV,k,lamda,t);
%对对角高频子带进行行与列的操作
dstCD=change(cD,k,lamda,t);
figure;
subplot(2,2,1);
imshow(cA);
title('低频图像');
subplot(2,2,2);
imshow(dstCH);
title('水平高频子带非线性扩散后的图像');
subplot(2,2,3);
imshow(dstCV);
title('垂直高频子带非线性扩散后的图像');
subplot(2,2,4);
imshow(dstCD);
title('对角高频子带非线性扩散后的图像');

3.1对图像的列进行两个方向的非线性扩散

function[dstI]=columchange(srcI,k,lamda,t)
%% 函数描述：对图像的列进行两个方向的非线性扩散
%% 输出参数
% dstI:非线性扩散后的图
%% 输入参数
% srcI:非线性扩散之前的图
% k:权衡系数
% lamda:控制平滑
% t:迭代次数
%% 函数代码
[r,c]=size(srcI);
for num=1:t
    for i=1:r
        for j=1:c
            if(i>1 && i<r && j>1 && j<c)
                left2right=srcI(i,j-1)-srcI(i,j);
                right2left=srcI(i,j+1)-srcI(i,j);
                c1=exp(-(left2right^2)/(k^2));
                c2=exp((-right2left^2)/(k^2));
                srcI(i,j)=srcI(i,j)+lamda*(left2right*c1+right2left*c2);
            end
        end
    end
end
dstI=srcI;

3.2对图像的行进行两个方向的非线性扩散

function[dstI]=rowchange(srcI,k,lamda,t)
%% 函数描述：对图像的行进行两个方向的非线性扩散
%% 输出参数
% dstI:非线性扩散后的图
%% 输入参数
% srcI:非线性扩散之前的图
% k:权衡系数
% lamda:控制平滑
% t:迭代次数
%% 函数代码
[r,c]=size(srcI);
for num=1:t
    for i=1:r
        for j=1:c
            if(i>1 && i<r && j>1 && j<c)
                up2down=srcI(i-1,j)-srcI(i,j);
                down2up=srcI(i+1,j)-srcI(i,j);
                c1=exp(-(up2down^2)/(k^2));
                c2=exp((-down2up^2)/(k^2));
                srcI(i,j)=srcI(i,j)+lamda*(up2down*c1+down2up*c2);
            end
        end
    end
end
dstI=srcI;

3.3对图像的行与列进行两个方向的非线性扩散

function[dstI]=change(srcI,k,lamda,t)
%% 函数描述：对图像的行与列进行两个方向的非线性扩散
%% 输出参数
% dstI:非线性扩散后的图
%% 输入参数
% srcI:非线性扩散之前的图
% k:权衡系数
% lamda:控制平滑
% t:迭代次数
%% 函数代码
[r,c]=size(srcI);
for num=1:t
    for i=1:r
        for j=1:c
            if(i>1 && i<r && j>1 && j<c)
                up2down=srcI(i-1,j)-srcI(i,j);
                down2up=srcI(i+1,j)-srcI(i,j);
                left2right=srcI(i,j-1)-srcI(i,j);
                right2left=srcI(i,j+1)-srcI(i,j);
                c1=exp(-(up2down^2)/(k^2));
                c2=exp((-down2up^2)/(k^2));
                c3=exp(-(left2right^2)/(k^2));
                c4=exp((-right2left^2)/(k^2));
                srcI(i,j)=srcI(i,j)+lamda*(up2down*c1+down2up*c2+left2right*c3+right2left*c4);
            end
        end
    end
end
dstI=srcI;

 4.逆小波变换

将各个高频子带进行去噪后与低频带再次拼接

%进行逆小波变换
tempI=iswt2(cA,dstCH,dstCV,dstCD,'haar');
figure;
imshow(tempI);
title('逆小波变换后的图像');

5.多尺度双边滤波

%多尺度双边滤波
r=5;
globalVar=3;
localVar=0.045;
dstI=doublefilter(tempI,r,globalVar,localVar);
figure;
imshow(dstI);
title('多尺度双边滤波后的结果图');

5.1 双边滤波函数

function[dstI]=doublefilter(srcI,r,globalVar,localVar)
%% 双边滤波函数
%% 输出参数
% dstI:双边滤波后的图像
%% 输入参数
% srcI:原始图像
% r:滤波半径
% globalVar:全局方差
% localVar:局部方差
%% 函数代码
%获取原始图像的宽高
[m,n]=size(srcI);
%计算空间权重
[x,y]=meshgrid(-r:r);
spaceW=exp(-(x.^2+y.^2)/(2*globalVar^2));
%填充原始图像
srcI_pad=padarray(srcI,[r r],'symmetric');
dstI=zeros(m,n);
for i=r+1:m+r
    for j=r+1:n+r
        %计算灰度权重
        temp=srcI_pad(i-r:i+r,j-r:j+r);
        grayW=exp(-(temp-srcI(i-r,j-r)).^2/(2*localVar^2));
        %计算最终权重
        finalW=spaceW.*grayW;
        ret=temp.*finalW;
        dstI(i-r,j-r)=sum(ret(:))/sum(finalW(:));
    end
end

6.均方差 峰值信噪比评价去噪能力

[MSE,PSNR]=Evaluation(srcI,dstI);
display(MSE);%均方差
display(PSNR);%峰值信噪比

6.1 评价函数

function [MSE,PSNR] = Evaluation( srcI,dstI )
%% 函数描述：评估去噪图像的效果
%% 输出参数：
% MSE：均方差
% PSNR：峰值信噪比
%% 输入参数：
% srcI:原始图片
% dstI：去噪后的图片
%% 函数代码
D=dstI-srcI;
MSE=sum(D(:).*D(:))/numel(dstI);
PSNR=10*log10(255^2/MSE);

7. 处理结果