实现option上下移动_Perona-Malik方程（各向同性非线性扩散实现图像滤波）

1. Perona-Malik方程

之前所使用的线性扩散图像滤波会对整幅图像“一视同仁”，不管是边缘纹理还是其他结构都会等同于噪声一起滤除。这会使图像过于平滑而且丢失很多高频信息。如果我们使用非线性扩散的方法，即：传导系数会随着图像的局部特征改变而改变，例如在图像平滑区域传导系数自适应增大可以更好的滤除噪声，在图像的边缘纹理区域系数自适应减小可以更好的保存图像边缘。

首先，PM扩散方程为：它的来源是将梯度的模值信息融合到传输系数中，并带入公式（5）可得：

其中

为梯度的模值：

，

为边缘函数，或者边缘停止函数。常用的边缘函数比如：

这里K是选定的常数，作用是控制g的下降速率。下图为p=1的条件下K=10, 20的下降曲线。

由于在图像的边缘中梯度的模值

取得极大值，而优化公式整体是求极小值的，为了避免在优化公式中引入负号，因此要引入一个单调递减而且是非负的边缘函数

，它可以保证在

取得极大值时，

取得局部极小值。它最早在图像分割和边缘检测问题中被提出，由于该函数的性质是可以将方程的解向着梯度变大的方向移动，最终停止在图像的边缘上，因此也被称为边缘停止函数。

2. P-M方程行为分析

首先将P-M方程简化为一维的情况分析：

其中

的下角标表示对某个变量求导。

称之为

影响函数。将（11）式中的g(-)函数带入可得：

下面是

的函数曲线图：

当p=1时，

从最大值1单调递减到0，而且始终保持正值，因此这种情况下的（12）式总保持正向扩散。
当p=2时，
这个不等式说明在梯度较小的图像平滑区域，即：

时，导数大于0为正向扩散。在梯度较大的边缘纹理区域，导数小于0为反向扩散。反向扩散的含义是：杂质将从浓度低的地方流向浓度高的地方，这意味着图像边缘的锐化。

在后面二维P-M方程的分析中，可以得出结论：当选用合适的边缘函数g时，P-M方程可以自适应的实现图像去噪和边缘增强的效果。

3. P-M方程的病态性质

数学研究表示（1）式给出的初值的P-M方程极有可能是病态的。如果定义下面这个能量泛函：

其中

是一个非负函数而且在0处取得的值为0。那么依据变分原理，最小化该能量泛函对应的梯度下降流为：

如果将

改写为:

，那么（5）式就和（1）式完全一致。也就是说P-M方程可以看做是（4）式能量泛函的梯度下降流。其边缘函数g由（4）式中的

函数决定。再与式（3）做比较可得函数的

的导数就是前面所定义的

：

4. P-M方程的正则化

一般来说，病态问题都是一个方程对应着多个解，因此需要用正则化的方法使它变成一个适定问题（well-posed）。有研究者提出一种新的改进方案为：

其中

表示方差为

的Gaussian函数，（7）式称之为正则化的P-M方程，或者叫CLMC模型。

5.P-M方程求解代码

搬运自：

数字图像处理，基于PM和Catte模型各向异性扩散的C++实现_EbowTang的练习场-CSDN博客blog.csdn.net

公式分析参考：

各向异性扩散PM模型原理与C++实现_cyh706510441的专栏-CSDN博客blog.csdn.net

function diff_im = anisodiff2D(im, num_iter, delta_t, kappa, option)
%ANISODIFF2D 经典PM模型的各向异性扩散
%   DIFF_IM = ANISODIFF2D(IM, NUM_ITER, DELTA_T, KAPPA, OPTION)
%   该函数执行灰度图像上的各向异性扩散(经典PM模型)，被认为是一个二维的网络结构的8个相邻节点的扩散传导。
%
%       参数描述:
%               IM       - 灰度图 (MxN).
%               NUM_ITER - 迭代次数
%               DELTA_T  - 积分常数 (0 <= delta_t <= 1/7).通常情况下，由于数值稳定性，此参数设置为它的最大值。
%               KAPPA    - 控制传导的梯度模阈值。控制平滑。
%               OPTION   - 传导系数函数选择（Perona & Malik提出）:
%                          1 - c(x,y,t) = exp(-(nablaI/kappa).^2),
%                              privileges high-contrast edges over low-contrast ones.
%                          2 - c(x,y,t) = 1./(1 + (nablaI/kappa).^2),
%                              privileges wide regions over smaller ones.
%
%       输出描述:
%                DIFF_IM - 具有最大尺度空间参数的（扩散）图像。
%
%   使用例子：
%   -------------
%   s = phantom(512) + randn(512);
%   num_iter = 15;
%   delta_t = 1/7;%越大越平滑
%   kappa = 30;%越大越平滑
%   option = 2;
%   ad = anisodiff2D(s,num_iter,delta_t,kappa,option);
%   figure, subplot 121, imshow(s,[]), subplot 122, imshow(ad,[]);% 转换输入图像类型为double.
im = double(im);% PDE（偏微分方程）的初始条件。
diff_im = im;% 中心像素距离。
dx = 1;
dy = 1;
dd = sqrt(2);% 二维卷积掩模-8个方向上的梯度差分。
hN = [0 1 0; 0 -1 0; 0 0 0];
hS = [0 0 0; 0 -1 0; 0 1 0];
hE = [0 0 0; 0 -1 1; 0 0 0];
hW = [0 0 0; 1 -1 0; 0 0 0];
hNE = [0 0 1; 0 -1 0; 0 0 0];
hSE = [0 0 0; 0 -1 0; 0 0 1];
hSW = [0 0 0; 0 -1 0; 1 0 0];
hNW = [1 0 0; 0 -1 0; 0 0 0];% 各向异性扩散
for t = 1:num_iter% 8个方向梯度差分. [imfilter(.,.,'conv') 也可以使用 conv2(.,.,'same')]nablaN = imfilter(diff_im,hN,'conv');nablaS = imfilter(diff_im,hS,'conv');   nablaW = imfilter(diff_im,hW,'conv');nablaE = imfilter(diff_im,hE,'conv');   nablaNE = imfilter(diff_im,hNE,'conv');nablaSE = imfilter(diff_im,hSE,'conv');   nablaSW = imfilter(diff_im,hSW,'conv');nablaNW = imfilter(diff_im,hNW,'conv'); % 扩散函数if option == 1cN = exp(-(nablaN/kappa).^2);cS = exp(-(nablaS/kappa).^2);cW = exp(-(nablaW/kappa).^2);cE = exp(-(nablaE/kappa).^2);cNE = exp(-(nablaNE/kappa).^2);cSE = exp(-(nablaSE/kappa).^2);cSW = exp(-(nablaSW/kappa).^2);cNW = exp(-(nablaNW/kappa).^2);elseif option == 2cN = 1./(1 + (nablaN/kappa).^2);cS = 1./(1 + (nablaS/kappa).^2);cW = 1./(1 + (nablaW/kappa).^2);cE = 1./(1 + (nablaE/kappa).^2);cNE = 1./(1 + (nablaNE/kappa).^2);cSE = 1./(1 + (nablaSE/kappa).^2);cSW = 1./(1 + (nablaSW/kappa).^2);cNW = 1./(1 + (nablaNW/kappa).^2);end% 离散偏微分方程的解决方案diff_im = diff_im + ...delta_t*(...(1/(dy^2))*cN.*nablaN + (1/(dy^2))*cS.*nablaS + ...(1/(dx^2))*cW.*nablaW + (1/(dx^2))*cE.*nablaE + ...(1/(dd^2))*cNE.*nablaNE + (1/(dd^2))*cSE.*nablaSE + ...(1/(dd^2))*cSW.*nablaSW + (1/(dd^2))*cNW.*nablaNW );% 迭代的警告fprintf('rIteration %dn',t);
end