图像处理（十一）图像分割(3)泛函能量LevelSet、snake分割

转自：https://blog.csdn.net/hjimce/article/details/45586727

一、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)是一个均值滤波的卷积核，因此如果要对算法用快速均值滤波，算法可提高五六倍的速度。具体算法代码实现如下：

[c++] view plain copy
 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                 原创文章，版权所有，转载请保留本行信息
--------------------- 
作者：hjimce 
来源：CSDN 
原文：https://blog.csdn.net/hjimce/article/details/45586727 
版权声明：本文为博主原创文章，转载请附上博文链接！