計算機視覺中的變分方法-擴散（Diffusion）
最近在看一個計算機視覺中的變分方法系列的視訊，是德國慕尼黑工大出的，講課老師是LSD-SLAM的作者Daniel Cremers，老師講得很清楚，看了還是很有收穫的。我已經變成Cremers大神的腦殘粉了，有興趣看視訊的戳這裡Variational Methods for Computer Vision
Diffusion equation:
擴散是一種物理過程，是讓空間中的物質的濃度分佈u(x,t)u(x,t)更加均勻一些。這個過程可以用兩個基礎的等式來描述：
Fick′slawFick′slaw : 空間物質的濃度的差別導致在濃度的負梯度方向上會有流jj。 這個也很好理解，意思就是說濃度高出的物質會往濃度低處擴散：
j=−g∇u
j=−g∇u
其中， gg是擴散係數（diffusivity），表示擴散過程的快慢
continuitycontinuity equationequation :
∂u∂t=−divj
∂u∂t=−divj
這裡，divj=∇⋅j=∂jx∂x+∂jy∂ydivj=∇·j=∂jx∂x+∂jy∂y 稱為散度。關於散度，其實在高等數學中有過介紹，通俗來講，對於空間場中一點，如果該點散度大於00，則表示該點向外擴散物質（好比是該點是水龍頭，向外流水）；如果該點散度等於00， 那就是擴散保持平衡，進多少出多少；如果散度小於00，那麼就說明該點在吸收物質（就像黑洞一樣吸收空間場中該點附近的物質）。
關於散度的更多資料，可以參見知乎上這個回答在影象處理中，散度 div 具體的作用是什麼
由上面兩個基本的等式，聯合起來就得到了今天要講的擴散方程（DiffusionDiffusion equationequation）
∂u∂t=div（g∇u）
∂u∂t=div（g∇u）
Example: Linear Diffusion Equation
下面以一維線性擴散方程為例來說明。
對於線性情況，g=1g=1:
∂tu=∂2tu
∂tu=∂t2u
初始條件：
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)=(G2t)∗f(x)=∫G2t(x−x′)f(x′)dx′
其中，Gσ=12πσ√exp(−x22σ2)Gσ=12πσexp(−x22σ2)，是一個高斯核，σ=2t−−√σ=2t
可以看到，高斯濾波其實是擴散的一種特例。但是我們都知道，用高斯濾波器對一個影象進行平滑濾波，由於高斯濾波器的各向同性，會使影象的邊緣細節都變模糊，有時候這不是我們想要的結果。
Nonlinear and Anisotropic Diffusion
一般形式下的擴散方程：
∂tu=div（g∇u）
∂tu=div（g∇u）
對影象濾波時，要想保持影象的邊緣細節，可以在影象邊緣資訊強的地方少擴散一些，那麼怎麼做呢？
我們用梯度的模來作為檢測邊緣的運算元|∇u|=u2x+u2y−−−−−−√|∇u|=ux2+uy2，那麼在邊緣處|∇u||∇u|的值就會比較大 ，然後再這些地方讓擴散速率變小，可以構造這樣的 gg:
g（|∇u|）=11+|∇u|2/λ2−−−−−−−−−−−√
g（|∇u|）=11+|∇u|2/λ2
其中，λ>0λ>0,稱為對比引數，在|∇u|>>λ|∇u|>>λ的區域，擴散速度接近於00，不受擴散的影響，所以可以保持該區域的細節。
關於這部分的詳細資料，可以參考影象處理的經典論文11Scale-space and edge detection using anisotropic diffusion
有限差分的實現
上面講了各向異性的擴散方程，現在就來說明一下如何程式設計實現。這部分內參考的是 Weickert, J Anisotropic diffusion in image processing裡的內容。
非線性擴散方程：
∂tu=∂x（g|∇u|∂xu）+∂y（g|∇u|∂yu）
∂tu=∂x（g|∇u|∂xu）+∂y（g|∇u|∂yu）
用差分來代替微分：∂tu=ut+1ij−utijτ∂tu=uijt+1−uijtτ
非線性擴散方程右邊第一項就可以表示為：
∂x（g∂xu）=(（g∂xu）ti+1/2,j−（g∂xu）ti−1/2,j)∂x（g∂xu）=(（g∂xu）i+1/2,jt−（g∂xu）i−1/2,jt)
=(gti+1/2,j(uti+1,j−uti,j)−gti−1/2,j(uti,j−uti−1,j))=(gi+1/2,jt(ui+1,jt−ui,jt)−gi−1/2,jt(ui,jt−ui−1,jt))
其中，gti+1/2,j=gi+1,jgi,j−−−−−−−√gi+1/2,jt=gi+1,jgi,j，說明只要這兩個畫素點處有一個的擴散速率gg為00，那麼插值得到的gti+1/2,jgi+1/2,jt就會為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程式碼示例關於程式碼實現的這部分內容，可以進一步參考這裡。
使用示例：
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
---------------------
作者：蝸牛一步一步往上爬
來源：CSDN
原文：https://blog.csdn.net/yc461515457/article/details/50847526
版權宣告：本文為博主原創文章，轉載請附上博文連結！
左邊是原圖，右邊是Anisotropic Diffusion結果圖
