计算机视觉中的变分方法-扩散（Diffusion）

最近在看一个计算机视觉中的变分方法系列的视频，是德国慕尼黑工大出的，讲课老师是LSD-SLAM的作者Daniel Cremers，老师讲得很清楚，看了还是很有收获的。我已经变成Cremers大神的脑残粉了，有兴趣看视频的戳这里Variational Methods for Computer Vision

Diffusion equation:

扩散是一种物理过程，是让空间中的物质的浓度分布u(x,t)u(x,t)更加均匀一些。这个过程可以用两个基础的等式来描述：

Fick′slawFick's law : 空间物质的浓度的差别导致在浓度的负梯度方向上会有流jj。 这个也很好理解，意思就是说浓度高出的物质会往浓度低处扩散：

j=−g∇u j = -g \nabla u

其中， gg是扩散系数（diffusivity），表示扩散过程的快慢

continuitycontinuity equation equation :

∂u∂t=−divj{{\partial u} \over {\partial t}} = - div j

这里，divj=∇⋅j=∂jx∂x+∂jy∂ydiv j =\nabla · j ={{\partial j_x} \over {\partial x}} + {{\partial j_y} \over {\partial y}} 称为散度。关于散度，其实在高等数学中有过介绍，通俗来讲，对于空间场中一点，如果该点散度大于00，则表示该点向外扩散物质（好比是该点是水龙头，向外流水）；如果该点散度等于00， 那就是扩散保持平衡，进多少出多少；如果散度小于00，那么就说明该点在吸收物质（就像黑洞一样吸收空间场中该点附近的物质）。

关于散度的更多资料，可以参见知乎上这个回答在图像处理中，散度 div 具体的作用是什么

由上面两个基本的等式，联合起来就得到了今天要讲的扩散方程（DiffusionDiffusion equationequation）

∂u∂t=div（g∇u）{{\partial u} \over {\partial t}} = div （g \nabla u）

Example: Linear Diffusion Equation

下面以一维线性扩散方程为例来说明。

对于线性情况，g=1g=1:

∂tu=∂2tu{{\partial _t u} } = \partial ^2 _t u

初始条件：u(x,t=0)=f(x) u(x,t=0) = f(x)

这个方程有唯一解：

u(x,t)=(G2t√)∗f(x)=∫G2t√(x−x′)f(x′)dx′u(x,t)=(G_{ \sqrt{2 t} })*f(x)= \int { G_{\sqrt{2t}} (x-x\prime )f(x\prime)dx\prime}

其中，Gσ=12πσ√exp(−x22σ2)G_\sigma = {{1}\over {\sqrt{2 \pi \sigma}} }exp(-{{x^2}\over {2\sigma^2}})，是一个高斯核，σ=2t−−√\sigma = \sqrt{2t}

可以看到，高斯滤波其实是扩散的一种特例。但是我们都知道，用高斯滤波器对一个图像进行平滑滤波，由于高斯滤波器的各向同性，会使图像的边缘细节都变模糊，有时候这不是我们想要的结果。

Nonlinear and Anisotropic Diffusion

一般形式下的扩散方程：

∂tu=div（g∇u）{\partial _t u} = div （g \nabla u）

对图像滤波时，要想保持图像的边缘细节，可以在图像边缘信息强的地方少扩散一些，那么怎么做呢？

我们用梯度的模来作为检测边缘的算子|∇u|=u2x+u2y−−−−−−√| \nabla u| = \sqrt{u_x^2+u^2_y}，那么在边缘处|∇u|| \nabla u|的值就会比较大 ，然后再这些地方让扩散速率变小，可以构造这样的 gg:

g（|∇u|）=11+|∇u|2/λ2−−−−−−−−−−−√g（| \nabla u|）={{1}\over{\sqrt {1+| \nabla u|^2/{\lambda^2}}}}

其中，λ>0\lambda >0,称为对比参数，在|∇u|>>λ| \nabla u| >>\lambda的区域，扩散速度接近于00，不受扩散的影响，所以可以保持该区域的细节。

关于这部分的详细资料，可以参考图像处理的经典论文1^1Scale-space and edge detection using anisotropic diffusion

有限差分的实现

上面讲了各向异性的扩散方程，现在就来说明一下如何编程实现。这部分内参考的是 Weickert, J Anisotropic diffusion in image processing里的内容。

非线性扩散方程：

∂tu=∂x（g|∇u|∂xu）+∂y（g|∇u|∂yu）{\partial _t u} = \partial_x（g |\nabla u|\partial_x u）+\partial_y（g |\nabla u|\partial_y u）

用差分来代替微分：∂tu=ut+1ij−utijτ\partial _t u = {{u^{t+1}_{ij} - u^{t}_{ij} } \over{\tau}}

非线性扩散方程右边第一项就可以表示为：

∂x（g∂xu）=(（g∂xu）ti+1/2,j−（g∂xu）ti−1/2,j) \partial_x（g \partial_x u）=(（g \partial_x u）^t_{i+1/2,j}-（g \partial_x u）^t_{i-1/2,j})

=(gti+1/2,j(uti+1,j−uti,j)−gti−1/2,j(uti,j−uti−1,j)) =( g^t_{i+1/2,j} (u^t_{i+1,j}-u^t_{i,j})-g^t_{i-1/2,j} (u^t_{i,j}-u^t_{i-1,j}))

其中，gti+1/2,j=gi+1,jgi,j−−−−−−−√g^t_{i+1/2,j} = \sqrt{g_{i+1,j}g_{i,j}}，说明只要这两个像素点处有一个的扩散速率gg为00，那么插值得到的gti+1/2,jg^t_{i+1/2,j} 就会为00，而不是去这两者的平均值。

关于这段差分实现的公式部分，需要说明的是扩散方程中对x,yx,y是进行了二阶差分，注意在上面公式中，第一次对xx方向差分选择的两个点是(i,j)(i,j)旁边的两个点(i+1/2,j)和(i−1/2,j)(i+1/2,j)和(i-1/2,j)，然后又进行了一次差分，得到的结果中，用到的像素点位置只有(i,j),(i+1,j),(i−1,j)(i,j),(i+1,j),(i-1,j),这样还是在一个3x33x3的窗口操作的，如果按照以前的第一次差分是右边的(i+1,j)(i+1,j)减左边的(i−1,j)(i-1,j)，那么结果就会出现(i+2,j),(i−2,j)(i+2,j),(i-2,j)项，最后就是相当于用了5x55x5的窗口，大的窗口对于细节的保持是不利的。

Anisotropic Diffusion Matlab代码示例

从网上找到的代码 matlab code

function diff_im = anisodiff2D(im, num_iter, delta_t, kappa, option)
%ANISODIFF2D Conventional anisotropic diffusion
%   DIFF_IM = ANISODIFF2D(IM, NUM_ITER, DELTA_T, KAPPA, OPTION) perfoms
%   conventional anisotropic diffusion (Perona & Malik) upon a gray scale
%   image. A 2D network structure of 8 neighboring nodes is considered for
%   diffusion conduction.
%
%       ARGUMENT DESCRIPTION:
%               IM       - gray scale image (MxN).
%               NUM_ITER - number of iterations.
%               DELTA_T  - integration constant (0 <= delta_t <= 1/7).
%                          Usually, due to numerical stability this
%                          parameter is set to its maximum value.
%               KAPPA    - gradient modulus threshold that controls the conduction.
%               OPTION   - conduction coefficient functions proposed by 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.
%
%       OUTPUT DESCRIPTION:
%                DIFF_IM - (diffused) image with the largest scale-space parameter.
%
%   Example
%   -------------
%   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,[])
%
% See also anisodiff1D, anisodiff3D.

% References:
%   P. Perona and J. Malik.
%   Scale-Space and Edge Detection Using Anisotropic Diffusion.
%   IEEE Transactions on Pattern Analysis and Machine Intelligence,
%   12(7):629-639, July 1990.
%
%   G. Grieg, O. Kubler, R. Kikinis, and F. A. Jolesz.
%   Nonlinear Anisotropic Filtering of MRI Data.
%   IEEE Transactions on Medical Imaging,
%   11(2):221-232, June 1992.
%
%   MATLAB implementation based on Peter Kovesi's anisodiff(.):
%   P. D. Kovesi. MATLAB and Octave Functions for Computer Vision and Image Processing.
%   School of Computer Science & Software Engineering,
%   The University of Western Australia. Available from:
%   <http://www.csse.uwa.edu.au/~pk/research/matlabfns/>.
%
% Credits:
% Daniel Simoes Lopes
% ICIST
% Instituto Superior Tecnico - Universidade Tecnica de Lisboa
% danlopes (at) civil ist utl pt
% http://www.civil.ist.utl.pt/~danlopes %
% May 2007 original version.

% Convert input image to double.
im = double(im);

% PDE (partial differential equation) initial condition.
diff_im = im;

% Center pixel distances.
dx = 1;
dy = 1;
dd = sqrt(2);

% 2D convolution masks - finite differences.
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];

% Anisotropic diffusion.
for t = 1:num_iter

% Finite differences. [imfilter(.,.,'conv') can be replaced by 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');

% Diffusion function.
if option == 1
cN = 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 == 2
cN = 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

% Discrete PDE solution.
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 );

% Iteration warning.
fprintf('\rIteration %d\n',t);
end

关于代码实现的这部分内容，可以进一步参考这里。

使用示例：

左边是原图，右边是Anisotropic Diffusion结果图

参考资料：

Variational Methods for Computer Vision - Lecture 4 (Prof. Daniel Cremers)

Pietro Perona and Jitendra Malik (July 1990). “Scale-space and edge detection using anisotropic diffusion”. IEEE Transactions on Pattern Analysis and Machine Intelligence 12 (7): 629–639.