各项异性扩散（Anisotropic diffusion）

各向异性扩散，也叫做P–M扩散，在图像处理和计算机视觉中广泛用于保持图像细节特征的同时减少噪声。

定义

有灰度图像$I(x,y)$，其各向异性扩散方程如下

$${{\partial I} \over {\partial t}} = div(c(x,y,t)\nabla I) = \nabla c\nabla I + c(x,y,t)\Delta I$$
其中 $\Delta$是Laplacian算子 , $\nabla$是梯度算子， $div$是散度， $c(x,y,t)$是扩散系数。控制着扩散速率，通常选取的图像梯度函数，这样在扩散时保护到图像边缘信息。这些由Perona和Malik在90年代初发现，他们提出两种扩散系数方程，也就是有名的P-M方程：
$$c(||\nabla I||) = {e^{ - {{(||\nabla I||/K)}^2}}}$$
和
$$c(||\nabla I||) = {1 \over {1 + {{({{||\nabla I||} \over K})}^2}}}$$
常数项K用来控制对边缘的灵敏度，通常经验选取或者用图像噪声相关的函数来表示。

原理

令$M$表示一副图像的光滑程度，那么上面的扩散方程可以用梯度下降方程最小化如下能量函数$E:M \to R$来表示

$$E[I] = {1 \over 2}\int_\Omega {g(||\nabla I(x)|{|^2})dx}$$
这里 $g:R \to R$是一个实数函数，稍后看到它本质上就是扩散系数。
对两边求其梯度，得到下式子
$$\nabla {E_I} = - div(g'(||\nabla I(x)|{|^2})\nabla I)$$
然后我们就可以使用梯度下降方程去降低图像的梯度，即平滑图像。
$${{\partial I} \over {\partial t}} = - \nabla {E_I} = div(g'(||\nabla I(x)||)\nabla I)$$
令 $c=g'$我们就得到了各向异性扩散方程。
好了，原理就解释到这么多了。若是想了解更多，可以看文章后面给出的更多阅读部分。

代码

% ANISODIFF - Anisotropic diffusion.
%
% Usage:
%  diff = anisodiff(im, niter, kappa, lambda, option)
%
% Arguments:
%         im     - input image
%         niter  - number of iterations.
%         kappa  - conduction coefficient 20-100 ?
%         lambda - max value of .25 for stability
%         option - 1 Perona Malik diffusion equation No 1
%                  2 Perona Malik diffusion equation No 2
%
% Returns:
%         diff   - diffused image.
%
% kappa controls conduction as a function of gradient.  If kappa is low
% small intensity gradients are able to block conduction and hence diffusion
% across step edges.  A large value reduces the influence of intensity
% gradients on conduction.
%
% lambda controls speed of diffusion (you usually want it at a maximum of
% 0.25)
%
% Diffusion equation 1 favours high contrast edges over low contrast ones.
% Diffusion equation 2 favours wide regions over smaller ones.

% Reference: 
% P. Perona and J. Malik. 
% Scale-space and edge detection using ansotropic diffusion.
% IEEE Transactions on Pattern Analysis and Machine Intelligence, 
% 12(7):629-639, July 1990.
%
% Peter Kovesi  
% School of Computer Science & Software Engineering
% The University of Western Australia
% pk @ csse uwa edu au
% http://www.csse.uwa.edu.au
%
% June 2000  original version.       
% March 2002 corrected diffusion eqn No 2.

function diff = anisodiff(im, niter, kappa, lambda, option)

if ndims(im)==3
  error('Anisodiff only operates on 2D grey-scale images');
end

im = double(im);
[rows,cols] = size(im);
diff = im;

for i = 1:niter
%  fprintf('\rIteration %d',i);

  % Construct diffl which is the same as diff but
  % has an extra padding of zeros around it.
  diffl = zeros(rows+2, cols+2);
  diffl(2:rows+1, 2:cols+1) = diff;

  % North, South, East and West differences
  deltaN = diffl(1:rows,2:cols+1)   - diff;
  deltaS = diffl(3:rows+2,2:cols+1) - diff;
  deltaE = diffl(2:rows+1,3:cols+2) - diff;
  deltaW = diffl(2:rows+1,1:cols)   - diff;

  % Conduction

  if option == 1
    cN = exp(-(deltaN/kappa).^2);
    cS = exp(-(deltaS/kappa).^2);
    cE = exp(-(deltaE/kappa).^2);
    cW = exp(-(deltaW/kappa).^2);
  elseif option == 2
    cN = 1./(1 + (deltaN/kappa).^2);
    cS = 1./(1 + (deltaS/kappa).^2);
    cE = 1./(1 + (deltaE/kappa).^2);
    cW = 1./(1 + (deltaW/kappa).^2);
  end

  diff = diff + lambda*(cN.*deltaN + cS.*deltaS + cE.*deltaE + cW.*deltaW);

%  Uncomment the following to see a progression of images
%  subplot(ceil(sqrt(niter)),ceil(sqrt(niter)), i)
%  imagesc(diff), colormap(gray), axis image

end
%fprintf('\n');

用途

感知边缘的滤波器在计算摄影学领域用途广泛，主要用于以下几个方面，仅举几例。
- 细节增强
- HDR色调映射
- 风格化
- 铅笔画
- 联合滤波
- 灰度图像彩色话等

效果

测试代码

 I = imread('lena.jpg');
 out = anisodiff(I,20,20,0.15,1);
 imshow(out/255);

参考阅读

WiKi百科英文
在图像处理中，散度 div 具体的作用是什么？
Scale-Space and Edge Detection Using Anisotropic Diffusion

转载请保留以下信息

作者	日期	联系方式
风吹夏天	2015年6月6日	wincoder@qq.com