影象處理（十一）影象分割(3)泛函能量LevelSet、snake分割
一、level set相關理論
基於水平集的影象分割演算法是一種進化版的Snake演算法，也是需要給定初始的輪廓曲線，然後根據泛函能量最小化，進行曲線演化。水平集的方法，用的是一種隱式函式的方法，這個演算法比較難理解，我一年前開始搞這個演算法的時候，雖然知道程式碼怎麼寫，但是它的原理推導完全不懂，因為這個演算法比較難理解，所以我這邊將講的稍微詳細一點。
跟傳統的snake演算法相比，思想完全不一樣，snake算法曲線演化的時候，是曲線上離散點顯示座標的位置更新移動，只要懂得能量最小化的曲線演化規則，就可以很快理解演算法，並寫出程式碼。然而水平集的方法，更新的不是曲線離散點的座標，而是更新整張圖片畫素點到曲線的有向距離場。因此演算法最關鍵的是理解這個距離場的更新規則，當然這個更新規則跟能量最小化相關。
開始這個演算法之前，我們需要非常熟悉，顯式二維的曲線與隱式曲線(水平集函式)的相互轉換公式。給定初始的輪廓曲線C，我們怎麼把它轉換成水平集函式，這個是實現演算法的第一步。水平集函式的定義：
公式1：
也就是說，如果給你一條初始封閉輪廓曲線C，進行水平集影象分割，我們需要寫的第一個函式就是計算影象的每個畫素點p(x,y)到曲線的最短距離d，如果該畫素點p位於曲線C的內部，那麼有向距離為-d；反之為d。這樣遍歷影象每個畫素點，每個畫素點都可以求得對應的有向距離u(x,y)。反過來，如果我已經知道了影象上每個畫素點的有向距離u(x,y)，那麼我要怎麼把這個隱式函式轉換成顯示函式呢？
其實很簡單，只要求出滿足u(x,y)=0 的畫素點，就是曲線上的點，因為如果該畫素點到曲線C的最短距離為0，那麼這個畫素點肯定在這條曲線上，據此我們就可以把所有滿足u(x,y)=0的畫素點全部提取出來，獲得這些畫素點的座標p(x,y)，而這些點便是曲線C的離散點，這樣就完成了從隱式距離場到顯式離散曲線的轉換。
據此我們可以得到演算法的大體流程：
輸入：給定離散的初始輪廓曲線C，待分割影象T
輸出：分割結果曲線C’
Algorithm:
Begin：
1、根據公式1，計算每個畫素點到離散曲線C的最短有向距離u(x,y)
2、根據影象梯度等資訊，對u(x,y)進行演化，使得其沿著能量最小化的方向演化，這個過程說的簡單一點，就是更新每個畫素點的u(x,y)值。
3、根據第2步的演化結果，遍歷每個畫素點(x,y)，判斷其水平集函式值是否為零。
If u(x,y)==0:
儲存畫素點座標(x,y)(因為這個點就是曲線C’上的點)
得到所有u(x,y)==0的點，就是最後我們想要的影象分割結果曲線C’
End
二、演算法實現：
這裡我選擇Mean Separation (MS) Energy能量最小化為例，講解區域性活動輪廓影象分割，具體的參考文獻為：《Localizing Region-Based Active Contours》。這裡為直接把文中最後演算法實現需要使用的公式，截出來，以便學習。
水平集總演化公式為：
其中：
dt是一個較小的數，選擇範圍為0.1~40都可以，當然迭代步長還是選擇的越小效果越好，就是需要的迭代次數越多。然後下面的公式為公式(17)對應的引數求解公式：
其中：
公式3是一個卷積核，上面的大部分過程的計算都涉及到用B(x,y)進行卷積，卷積半徑在我的專案使用的時候，我是選擇R=8，而且B(x,y)是一個均值濾波的卷積核，因此如果要對演算法用快速均值濾波，演算法可提高五六倍的速度。具體演算法程式碼實現如下：
function seg = local_AC_MS(Img,mask_init,rad,alpha,num_it,epsilon)
%引數Img為灰度影象
%mask_init為初始曲線的一個mask，曲線上的點座標在mask中的值為1
%rad 為卷積半徑
%alpha為公式6中的λ
%num_it為迭代次數
%epsilon為公式1中的ε
phi0 = mask2phi(mask_init);%從顯式曲線轉換成水平集函式
phi = phi0;%計算距離場，即計算所有的畫素點位置（x,y）到曲線的最短距離，這就是所謂的水平集函式
B0 = ones(2*rad+1,2*rad+1); %mask，卷積核
KI=conv2(Img,B0,'same'); %對影象Img用卷積核B0進行卷積
KONE=conv2(ones(size(Img)),B0,'same');
%下面計算的是文獻的公式17，即使用公式17對曲線進行演化，所有的計算都是為了計算公式17
for ii = 1:num_it%開始迭代
mask = Heaviside2(phi,epsilon);%計算文獻公式1
I=Img.*mask;
temp1=conv2(mask,B0,'same'); %文獻的公式18
temp2=conv2(I,B0,'same');
c1=temp2./(temp1);
c2=(KI-temp2)./(KONE-temp1);
A1 = temp1;%文獻的公式18
A2 = conv2(1-mask,B0,'same');%文獻公式19
D = (A1.*A2+eps);
term1 = (A2-A1)./D;
term2 = (A2.*c1.^2-A1.*c2.^2)./D;
term3 = (A2.*c1-A1.*c2)./D;
%計算總公式的前半部分
dataForce = conv2(term1.*Dirac2(phi,epsilon),B0,'same').*Img.*Img + conv2(term2.*Dirac2(phi,epsilon),B0,'same')-2.*Img.*conv2(term3.*Dirac2(phi,epsilon),B0,'same'); %%% During the implementation, Img should be separated out of the filtering operation!!!
dataForce = dataForce/max(abs(dataForce(:)));
%計算總公式的後半部分，即水平集的散度
curvature = curvature_central(phi);
%總公式
dphi = Dirac2(phi,epsilon).*(-dataForce + alpha*curvature);
dt = 1/(max(abs(dphi(:)))+eps);%時間步長，該引數可人為設定為恆定引數
%曲線演化公式，即完成曲線的迭代
phi = phi + dt.*dphi;
%繪製曲線，（x，y）的值為0的點即為曲線上的點
if(mod(ii,10) == 0)
showCurveAndPhi(Img,phi,ii); %繪製值為0的等值線
end
end
seg = (phi>=0);
%為了保證水平集的光滑性，需要對水平集進行重新計算，保證水平集的梯度模為1
%具體計算公式見
function D = sussman(D, dt)
%前後差分
a = D - shiftR(D); % 後差分
b = shiftL(D) - D; % 前差分
c = D - shiftD(D);
d = shiftU(D) - D;
a_p = a; a_n = a; % a+ and a-
b_p = b; b_n = b;
c_p = c; c_n = c;
d_p = d; d_n = d;
a_p(a < 0) = 0;
a_n(a > 0) = 0;
b_p(b < 0) = 0;
b_n(b > 0) = 0;
c_p(c < 0) = 0;
c_n(c > 0) = 0;
d_p(d < 0) = 0;
d_n(d > 0) = 0;
dD = zeros(size(D));
D_neg_ind = find(D < 0);
D_pos_ind = find(D > 0);
dD(D_pos_ind) = sqrt(max(a_p(D_pos_ind).^2, b_n(D_pos_ind).^2) ...
+ max(c_p(D_pos_ind).^2, d_n(D_pos_ind).^2)) - 1;
dD(D_neg_ind) = sqrt(max(a_n(D_neg_ind).^2, b_p(D_neg_ind).^2) ...
+ max(c_n(D_neg_ind).^2, d_p(D_neg_ind).^2)) - 1;
D = D - dt .* sussman_sign(D) .* dD;
function shift = shiftD(M)
shift = shiftR(M')';
function shift = shiftL(M)
shift = [ M(:,2:size(M,2)) M(:,size(M,2)) ];
function shift = shiftR(M)
shift = [ M(:,1) M(:,1:size(M,2)-1) ];
function shift = shiftU(M)
shift = shiftL(M')';
function S = sussman_sign(D)
S = D ./ sqrt(D.^2 + 1);
%計算有向距離場函式
function phi = mask2phi(init_a)
phi=bwdist(init_a)-bwdist(1-init_a)+im2double(init_a);
phi = -double(phi);
%水平集提取函式，也就是把隱式函式轉換成顯式函式，所得簡單一點，就是提取值為0的等高線
function showCurveAndPhi(I, phi, i)
imshow(I,'initialmagnification',200,'displayrange',[ ]);
hold on; contour(phi, [0 0], 'y','LineWidth',2);
hold off; title([num2str(i) ' Iterations']); drawnow;
%文獻公式2
function f = Dirac2(x, sigma)
f = (sigma/pi)./(sigma^2+x.^2);
%計算文獻公式1
function f = Heaviside2(x, epsilon)
f = 0.5*(1+2/pi*atan(x./epsilon));
%散度計算
function k = curvature_central(u)
[ux,uy] = gradient(u);
normDu = sqrt(ux.^2+uy.^2+1e-10);
Nx = ux./normDu;
Ny = uy./normDu;
[nxx,junk] = gradient(Nx);
[junk,nyy] = gradient(Ny);
k = nxx+nyy; 個人觀點：我覺得傳統的snake演算法，只能用爛來解釋，基本上以遇到一點噪聲，就不行了，分割精度真不是一般的差。而水平集的分割只能用高大上來形容，可以進行自動分裂合併等，就是速度有點慢，因為每次都要對水平集函式進行更新，更新一次就相當於遍歷一張圖片，因此速度可想而知，當然還有很多加速版的水平集方法，有待測試學習。本文地址：http://blog.csdn.net/hjimce/article/details/45586727 作者：hjimce 聯絡qq：1393852684 更多資源請關注我的部落格：http://blog.csdn.net/hjimce 原創文章，版權所有，轉載請保留本行資訊