最近在看一個計算機視覺中的變分方法系列的視訊,是德國慕尼黑工大出的,講課老師是LSD-SLAM的作者Daniel Cremers,老師講得很清楚,看了還是很有收獲的。我已經變成Cremers大神的腦殘粉了,有興趣看視訊的戳這裡Variational Methods for Computer Vision
Diffusion equation:
擴散是一種實體過程,是讓空間中的物質的濃度分布 u(x,t) 更加均勻一些。這個過程可以用兩個基礎的等式來描述:
-
Fick′slaw : 空間物質的濃度的差别導緻在濃度的負梯度方向上會有流 j 。 這個也很好了解,意思就是說濃度高出的物質會往濃度低處擴散:
j=−g∇u
其中, g 是擴散系數(diffusivity),表示擴散過程的快慢
-
continuity equation :
∂u∂t=−divj
這裡, divj=∇⋅j=∂jx∂x+∂jy∂y 稱為散度。關于散度,其實在高等數學中有過介紹,通俗來講,對于空間場中一點,如果該點散度大于 0 ,則表示該點向外擴散物質(好比是該點是水龍頭,向外流水);如果該點散度等于0, 那就是擴散保持平衡,進多少出多少;如果散度小于 0 ,那麼就說明該點在吸收物質(就像黑洞一樣吸收空間場中該點附近的物質)。
關于散度的更多資料,可以參見知乎上這個回答在圖像進行中,散度 div 具體的作用是什麼
由上面兩個基本的等式,聯合起來就得到了今天要講的擴散方程(Diffusion equation )
∂u∂t=div(g∇u)
Example: Linear Diffusion Equation
下面以一維線性擴散方程為例來說明。
對于線性情況, g=1 :
∂tu=∂2tu
初始條件: u(x,t=0)=f(x)
這個方程有唯一解:
u(x,t)=(G2t√)∗f(x)=∫G2t√(x−x′)f(x′)dx′
其中, Gσ=12πσ√exp(−x22σ2) ,是一個高斯核, σ=2t−−√
可以看到,高斯濾波其實是擴散的一種特例。但是我們都知道,用高斯濾波器對一個圖像進行平滑濾波,由于高斯濾波器的各向同性,會使圖像的邊緣細節都變模糊,有時候這不是我們想要的結果。
Nonlinear and Anisotropic Diffusion
一般形式下的擴散方程:
∂tu=div(g∇u)
對圖像濾波時,要想保持圖像的邊緣細節,可以在圖像邊緣資訊強的地方少擴散一些,那麼怎麼做呢?
我們用梯度的模來作為檢測邊緣的算子 |∇u|=u2x+u2y−−−−−−√ ,那麼在邊緣處 |∇u| 的值就會比較大 ,然後再這些地方讓擴散速率變小,可以構造這樣的 g :
g(|∇u|)=11+|∇u|2/λ2−−−−−−−−−−−√
其中, λ>0 ,稱為對比參數,在 |∇u|>>λ 的區域,擴散速度接近于 0 ,不受擴散的影響,是以可以保持該區域的細節。
關于這部分的詳細資料,可以參考圖像處理的經典論文1Scale-space and edge detection using anisotropic diffusion
有限差分的實作
上面講了各向異性的擴散方程,現在就來說明一下如何程式設計實作。這部分内參考的是 Weickert, J Anisotropic diffusion in image processing裡的内容。
非線性擴散方程:
∂tu=∂x(g|∇u|∂xu)+∂y(g|∇u|∂yu)
用差分來代替微分: ∂tu=ut+1ij−utijτ
非線性擴散方程右邊第一項就可以表示為:
∂x(g∂xu)=((g∂xu)ti+1/2,j−(g∂xu)ti−1/2,j)
=(gti+1/2,j(uti+1,j−uti,j)−gti−1/2,j(uti,j−uti−1,j))
其中, gti+1/2,j=gi+1,jgi,j−−−−−−−√ ,說明隻要這兩個像素點處有一個的擴散速率 g 為0,那麼插值得到的 gti+1/2,j 就會為 0 ,而不是去這兩者的平均值。
關于這段差分實作的公式部分,需要說明的是擴散方程中對x,y是進行了二階差分,注意在上面公式中,第一次對 x 方向差分選擇的兩個點是(i,j)旁邊的兩個點 (i+1/2,j)和(i−1/2,j) ,然後又進行了一次差分,得到的結果中,用到的像素點位置隻有 (i,j),(i+1,j),(i−1,j) ,這樣還是在一個 3x3 的視窗操作的,如果按照以前的第一次差分是右邊的 (i+1,j) 減左邊的 (i−1,j) ,那麼結果就會出現 (i+2,j),(i−2,j) 項,最後就是相當于用了 5x5 的視窗,大的視窗對于細節的保持是不利的。Anisotropic Diffusion Matlab代碼示例
從網上找到的代碼 matlab codefunction 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.