图像处理（十一）图像分割(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                 原创文章，版权所有，转载请保留本行信息