主要是用来平滑图像的，克服了高斯模糊的缺陷，各向异性扩散在平滑图像时是保留图像边缘的(和双边滤波很像)。

通常我们有将图像看作矩阵的，看作图的，看作随机过程的，记得过去还有看作力场的。

这次新鲜，将图像看作热量场了。每个像素看作热流，根据当前像素和周围像素的关系，来确定是否要向周围扩散。比如某个邻域像素和当前像素差别较大，则代表这个邻域像素很可能是个边界，那么当前像素就不向这个方向扩散了，这个边界也就得到保留了。

先看下效果吧：

具体的推导公式都是热学上的，自己也不太熟悉，感兴趣的可以去看原论文，引用量超7000吶。

我这里只介绍一下最终结论用到的公式。

主要迭代方程如下：

I就是图像了，因为是个迭代公式，所以有迭代次数t。

四个散度公式是在四个方向上对当前像素求偏导，news就是东南西北嘛，公式如下：

而cN/cS/cE/cW则代表四个方向上的导热系数，边界的导热系数都是小的。公式如下：

最后整个公式需要先前设置的参数主要有三个，迭代次数t，根据情况设置；导热系数相关的k，取值越大越平滑，越不易保留边缘；lambda同样也是取值越大越平滑。

最后是matlab代码：

clear all;

close all;

clc;

k=15; %导热系数,控制平滑

lambda=0.15; %控制平滑

N=20; %迭代次数

img=double(imread('lena.jpg'));

imshow(img,[]);

[m n]=size(img);

imgn=zeros(m,n);

for i=1:N

for p=2:m-1

for q=2:n-1

%当前像素的散度，对四个方向分别求偏导，局部不同方向上的变化量，

%如果变化较多，就证明是边界，想方法保留边界

NI=img(p-1,q)-img(p,q);

SI=img(p+1,q)-img(p,q);

EI=img(p,q-1)-img(p,q);

WI=img(p,q+1)-img(p,q);

%四个方向上的导热系数，该方向变化越大，求得的值越小，从而达到保留边界的目的

cN=exp(-NI^2/(k*k));

cS=exp(-SI^2/(k*k));

cE=exp(-EI^2/(k*k));

cW=exp(-WI^2/(k*k));

imgn(p,q)=img(p,q)+lambda*(cN*NI+cS*SI+cE*EI+cW*WI); %扩散后的新值

end

end

img=imgn; %整个图像扩散完毕，用已扩散图像的重新扩散。

end

figure;

imshow(imgn,[]);

VC代码(借助于opencv)：

int k = 15;

float lambda = 0.25;

int N = 20;

int m = im_eye.rows;

int n = im_eye.cols;

Mat imgn(m, n, CV_32FC1, Scalar(0));

for (int i = 0; i < N; i++)

{

for (int p = 1; p < m-1; p++)

{

for (int q = 1; q < n-1; q++)

{

float NI = diffimg.at(p - 1, q) - diffimg.at(p, q);

float SI = diffimg.at(p + 1, q) - diffimg.at(p, q);

float EI = diffimg.at(p, q - 1) - diffimg.at(p, q);

floatWI = diffimg.at(p, q + 1) - diffimg.at(p, q);

floatcN = exp(-pow(NI ,2) / (k*k));

float cS = exp(-pow(SI, 2) / (k*k));

floatcE = exp(-pow(EI, 2) / (k*k));

floatcW = exp(-pow(WI, 2) / (k*k));

imgn.at(p, q) = diffimg.at(p, q) + lambda*(cN*NI + cS*SI + cE*EI + cW*WI);

}

}

diffimg = imgn;

}