图像滤波---各项异性扩散滤波使用笔记及代码

一、文章内容介绍

1、文章首先介绍各项异性扩散滤波以及各项同性滤波
2、接着简要的说明各项异性滤波算法的原理和功能
3、然后对各项异性滤波算法的几个参数的作用及调试进行详解
4、最后附上各项异性滤波matlab的代码

二、各项异性扩散滤波和各项同性滤波

1、各项同性滤波

在做图像处理时，我们接触到很多滤波算法，耳熟能详的有高斯滤波、均值滤波以及中值滤波。这些滤波
算法中，对于每一个像素点，算法都一视同仁，或使用相同的系数，或使用相同的邻域规则来对像素值进行处理。这些算法不考虑该像素的梯度方向，一视同仁，处理效果在各个方向上是一样的，对于整幅图像的处理是均匀的，图像的平滑不会因梯度而有所偏向，这样的滤波叫各项同性滤波（Isotropic Filtering）。

2、各项异性扩散滤波

与各项同性扩散滤波相对的即是各项异性扩散滤波（Anisotropic Filtering），各项异性扩散滤波是一种保边滤波器，在不同的梯度方向上，像素值做不同程度的改变，区分出图像的边缘信息，在去除噪声的同时保留边缘。
各项异性扩散叫P-M扩散，论文原文：Scale-space and edge detection using anisotropic diffusion。

3、各项异性和各项同性的对比

各项同性滤波处理是均匀的，不考虑方向，适合于简单的去噪和模糊处理。
各项异性滤波处理是方向依赖的，能够更好地保留边缘和细节，适合于复杂的图像去噪，尤其是在需要保留边缘信息的情况下。

三、各项异性扩散滤波的原理介绍

各项异性扩散滤波的本质是根据图像的梯度信息来调控扩散系数。
各项异性指物体的全部或者部分物理、化学性质随方向的不同而有所变化的特性。对于整张图像，将其看做一个热量场，将每个像素看做一个热流，对于边缘像素，我们不想让其平滑掉，所以边缘的邻域像素，热流的扩散要很弱，除此之外，其他的地方，顺着梯度的方向进行扩散，扩散过的图像区域就会变得平滑，如此在保留边缘的同时，平滑了噪声区域。

四、各项异性扩散滤波公式

滤波公式如下：
$$I_{t+1}=I_t+\lambda \ast (c{N}_{xy}\ast D_N(I_t)+c{S}_{xy}\ast D_S(I_t)+c{E}_{xy}\ast D_E(I_t)+c{W}_{xy}\ast D_W(I_t))$$
其中I表示图像，x,y表示像素位置，t代表迭代次数，$\lambda$表示扩散速率。$D_N$、$D_S$、$D_E$、$D_W$分别表示南、北、西、东，即上下左右四个方向的梯度，梯度计算如下：
$$D_N(I_t)=I_t(x,y-1)-I_t(x,y)$$$$D_S(I_t)=I_t(x,y+1)-I_t(x,y)$$$$D_E(I_t)=I_t(x-1,y)-I_t(x,y)$$$$D_W(I_t)=I_t(x+1,y)-I_t(x,y)$$
Nxy、Sxy、Exy、Wxy分别表示对应方向上的扩散系数，扩散系数公式如下：
$$c{N}_{xy}=e^{-\frac{{D_N(I_t)}^2}{K^2}}$$$$c{S}_{xy}=e^{-\frac{{D_S(I_t)}^2}{K^2}}$$$$c{E}_{xy}=e^{-\frac{{D_E(I_t)}^2}{K^2}}$$$$c{W}_{xy}=e^{-\frac{{D_W(I_t)}^2}{K^2}}$$
其中$K$表示对应方向的扩散灵敏度。

五、公式中的参数使用说明

1、扩散速率$\lambda$

（1）参数说明

从公式中显而易见，$\lambda$控制每次迭代中像素值的变化量。
理论上，$\lambda$较大每个像素的扩散速度会比较快，因为每次像素值改变的增量会比较大，会导致图像的细节丢失。

（2）效果展示

在N(迭代次数)和K固定的情况下，$\lambda$从0.1~0.9的去噪效果图，0.1时细节保留最完整，0.3时细节丢失严重，0.5 ~0.7细节丢失程度依次增加。

（3）调节指南

从效果图分析，调解时，固定另外两个变量，从0.1开始逐步增加或减少，步长可设置为0.05、0.1、0.2，观察图像的去噪效果和细节保留情况，找到最佳的平衡点，比如，上述图像的滤波，平衡点应在0.1 ~0.3之间。

2、扩散敏感度$K$

（1）参数说明

从公式看，$K$是指数函数$e$的一个参数，函数形式为$y=e^{-\frac{1}{K^2}\ast x^2}$，当$\frac{1}{K^2}=0.5$时，函数图像如下：

当$\frac{1}{K^2}=2$时，函数图像如下：

从两幅图像来看，$K$越大，曲线变换越平缓，表明系数对梯度的响应越不敏感，各项异性的程度就越低，图像应该就越模糊。

（2）效果展示

固定$\lambda =0.12$，N=20时，$K$从10增加到40，从图上看，$K$越大，平滑的越厉害，细节丢失越严重。

（3）调节指南

对于$K$，通常可以选择从10开始调节，增加或减少，每次步长5、10、15、20尝试。
如果图像中细节多噪声少，K的值应该小一些，否则就大一些。

六、完整MATLAB代码

运行环境为Windows10和Matlab2023a
run.m

clc 
clear
%读取图像
img = imread('lisaGrayNoise.bmp');
%设置系数
lambda = 0.12;
K = 40;
N = 20;
%执行各项异性扩散滤波
imgADF = ImageADF(img,lambda,K,N);
figure;
imshow(img);
title('原图');
figure;
imshow(imgADF);
title('各向异性扩散滤波结果K=40');

ImageADF.m

function imgADF = ImageADF(img,lambda,K,N)
[m, n] = size(img);
imgADF = zeros(m, n);
img = double(img);
for t = 1 : N
    for i = 2 : m-1
        for j = 2 : n-1
            NI = img(i-1, j) - img(i, j);
            SI = img(i+1, j) - img(i, j);
            EI = img(i, j-1) - img(i, j);
            WI = img(i, j+1) - img(i, j);
            cN = exp(-NI^2/K^2);
            cS = exp(-SI^2/K^2);
            cE = exp(-EI^2/K^2);
            cW = exp(-WI^2/K^2);
            imgADF(i, j) = img(i, j) + lambda*(NI*cN + SI*cS + EI*cE + WI*cW);
        end
    end
    img = imgADF;
end
imgADF = img;
imgADF = uint8(imgADF);
end