各项异性扩散滤波

1990年Perona和Malik提出了各向异性扩散方程( 原论文)根据不同方向上的梯度来确定扩散系数,在抑制噪声的同时提高了保留细节的能力

1.关键词

各向异性扩散滤波（Anisotropic Filter，也叫P–M扩散），同中值滤波、高斯滤波一样，也是一种滤除图像噪音的方法。可以用于边缘检测或提取之前的预处理，去除无关噪声，增加边缘的有效性。某些基于各向异性滤波的边缘提取算法，在滤波的同时还能兼顾边缘保持。

各向异性，指物体的全部或者部分物理、化学性质随方向的不同而有所变化的特性。例如石墨单晶的电导率在不同方向的差异可达数千倍，又如天文学上，宇宙微波背景辐射亦拥有些微的各向异性。许多的物理量都具有各向异性，如弹性模量、电导率、在酸中的溶解速度等。

各向异性扩散，是有选择性的平滑（扩散）过程，这种选择是指方向选择的非均匀性，这种平滑过程在均匀的区域不受限制，而在跨越边界部分被抑制，因此在平滑噪声的同时保持图像的边缘特征。

转载于：https://blog.csdn.net/qq_32391345/article/details/106750055

回忆一下散度的概念：通量就是流过一个固定表面的流量（可以想象成水管截面），散度就是单位面积的通量。

2.定义

有灰度图像I(x,y)I(x,y)，其各向异性扩散方程如下：
$$\frac{\partial I}{\partial t}=\mathrm{div}(c(x,y,t)\nabla I)=\nabla c\nabla I+c(x,y,t)\Delta I$$

其中 Δ是Laplacian算子 , ∇是梯度算子， div是散度， c(x,y,t)是扩散系数。控制着扩散速率，通常选取的图像梯度函数，这样在扩散时保护到图像边缘信息。这些由Perona和Malik在90年代初发现，他们提出两种扩散系数方程，也就是有名的P-M方程：
$$c(\Vert \nabla I\Vert )=e^{-(\Vert \nabla I\Vert /K{)}^2}$$

和
$$c(\Vert \nabla I\Vert )=\frac{1}{1+{\left(\frac{\Vert \nabla I\Vert }{K}\right)}^2}$$
常数项K用来控制对边缘的灵敏度，通常经验选取或者用图像噪声相关的函数来表示。

3. 原理

转载于

4.滤波公式

$$I_{t+1}=I_t+\lambda \left(c{N}_{x,y}{\nabla }_N\left(I_t\right)+c{S}_{x,y}{\nabla }_S\left(I_t\right)+c{E}_{x,y}{\nabla }_E\left(I_t\right)+c{W}_{x,y}{\nabla }_W\left(I_t\right)\right)$$
I就是图像，因为是个迭代公式，所以有迭代次数t。

四个散度公式是在四个方向上对当前像素求偏导，NSEW就是东南西北，公式如下：
$$\begin{array}{rl}& {\nabla }_N\left(I_{x,y}\right)=I_{x,y-1}-I_{x,y}\\ & {\nabla }_S\left(I_{x,y}\right)=I_{x,y+1}-I_{x,y}\\ & {\nabla }_E\left(I_{x,y}\right)=I_{x-1,y}-I_{x,y}\\ & {\nabla }_W\left(I_{xy}\right)=I_{x+1,y}-I_{x,y}\end{array}$$
而cN/cS/cE/cW则代表四个方向上的导热系数，边界的导热系数都是小的。公式如下：$$\begin{array}{rl}c{N}_{x,y}& =\exp \left(-{\Vert {\nabla }_N(I)\Vert }^2/k^2\right)\\ c{S}_{x,y}& =\exp \left(-{\Vert {\nabla }_s(I)\Vert }^2/k^2\right)\\ c{E}_{x,y}& =\exp \left(-{\Vert {\nabla }_E(I)\Vert }^2/k^2\right)\\ c{W}_{x,y}& =\exp \left(-{\Vert {\nabla }_W(I)\Vert }^2/k^2\right)\end{array}$$

5.代码

function newPic = PM(pic)
[m n] = size(pic);
lambda = 0.20;
K = 24;
N = 40;
newPic = zeros(m, n);
pic = double(pic);
for t = 1 : N
    for i = 2 : m-1
        for j = 2 : n-1
            NI = pic(i-1, j) - pic(i, j);
            SI = pic(i+1, j) - pic(i, j);
            EI = pic(i, j-1) - pic(i, j);
            WI = pic(i, j+1) - pic(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);
            %fprintf('%d %d\n', i, j);
            newPic(i, j) = pic(i, j) + lambda*(NI*cN + SI*cS + EI*cE + WI*cW);
        end;
    end
    pic = newPic;
end;
newPic = pic;
newPic = uint8(newPic);
end

clc 
clear
%%各项异性扩散
pic = imread('2.jpg');
pic = rgb2gray(pic);
picNew = PM(pic);
figure;
imshow(pic);
title('原图');
figure;
imshow(picNew);
title('各向异性扩散');
% imwrite(picNew,'3.jpg');