计算机视觉
图像处理

图像处理(九)图像分割(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)是一个均值滤波的卷积核,因此如果要对算法用快速均值滤波,算法可提高五六倍的速度。具体算法代码实现如下:

  1.  function seg = local_AC_MS(Img,mask_init,rad,alpha,num_it,epsilon)
  2. %参数Img为灰度图像
  3. %mask_init为初始曲线的一个mask,曲线上的点坐标在mask中的值为1
  4. %rad 为卷积半径
  5. %alpha为公式6中的λ
  6. %num_it为迭代次数
  7. %epsilon为公式1中的ε
  8. phi0 = mask2phi(mask_init);%从显式曲线转换成水平集函数
  9. phi = phi0;%计算距离场,即计算所有的像素点位置(x,y)到曲线的最短距离,这就是所谓的水平集函数
  10. B0 = ones(2*rad+1,2*rad+1);  %mask,卷积核
  11. KI=conv2(Img,B0,‘same’);  %对图像Img用卷积核B0进行卷积
  12. KONE=conv2(ones(size(Img)),B0,‘same’);
  13. %下面计算的是文献的公式17,即使用公式17对曲线进行演化,所有的计算都是为了计算公式17
  14. for ii = 1:num_it%开始迭代
  15. mask = Heaviside2(phi,epsilon);%计算文献公式1
  16. I=Img.*mask;
  17. temp1=conv2(mask,B0,‘same’);   %文献的公式18
  18. temp2=conv2(I,B0,‘same’);
  19. c1=temp2./(temp1);
  20. c2=(KI-temp2)./(KONE-temp1);
  21. A1 = temp1;%文献的公式18
  22. A2 = conv2(1-mask,B0,‘same’);%文献公式19
  23. D = (A1.*A2+eps);
  24. term1 = (A2-A1)./D;
  25. term2 = (A2.*c1.^2-A1.*c2.^2)./D;
  26. term3 = (A2.*c1-A1.*c2)./D;
  27. %计算总公式的前半部分
  28. 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!!!
  29. dataForce = dataForce/max(abs(dataForce(:)));
  30. %计算总公式的后半部分,即水平集的散度
  31. curvature = curvature_central(phi);
  32. %总公式
  33. dphi = Dirac2(phi,epsilon).*(-dataForce + alpha*curvature);
  34. dt = 1/(max(abs(dphi(:)))+eps);%时间步长,该参数可人为设置为恒定参数
  35. %曲线演化公式,即完成曲线的迭代
  36. phi = phi + dt.*dphi;
  37. %绘制曲线,(x,y)的值为0的点即为曲线上的点
  38.     if(mod(ii,10) == 0)
  39.       showCurveAndPhi(Img,phi,ii);  %绘制值为0的等值线
  40.     end
  41. end
  42. seg = (phi>=0);
  43. %为了保证水平集的光滑性,需要对水平集进行重新计算,保证水平集的梯度模为1
  44. %具体计算公式见
  45. function D = sussman(D, dt)
  46.   %前后差分
  47.   a = D – shiftR(D); % 后差分
  48.   b = shiftL(D) – D; % 前差分
  49.   c = D – shiftD(D);
  50.   d = shiftU(D) – D;
  51.   a_p = a;  a_n = a; % a+ and a-
  52.   b_p = b;  b_n = b;
  53.   c_p = c;  c_n = c;
  54.   d_p = d;  d_n = d;
  55.   a_p(a < 0) = 0;
  56.   a_n(a > 0) = 0;
  57.   b_p(b < 0) = 0;
  58.   b_n(b > 0) = 0;
  59.   c_p(c < 0) = 0;
  60.   c_n(c > 0) = 0;
  61.   d_p(d < 0) = 0;
  62.   d_n(d > 0) = 0;
  63.   dD = zeros(size(D));
  64.   D_neg_ind = find(D < 0);
  65.   D_pos_ind = find(D > 0);
  66.   dD(D_pos_ind) = sqrt(max(a_p(D_pos_ind).^2, b_n(D_pos_ind).^2) …
  67.                        + max(c_p(D_pos_ind).^2, d_n(D_pos_ind).^2)) – 1;
  68.   dD(D_neg_ind) = sqrt(max(a_n(D_neg_ind).^2, b_p(D_neg_ind).^2) …
  69.                        + max(c_n(D_neg_ind).^2, d_p(D_neg_ind).^2)) – 1;
  70.   D = D – dt .* sussman_sign(D) .* dD;
  71. function shift = shiftD(M)
  72.   shift = shiftR(M‘)’;
  73. function shift = shiftL(M)
  74.   shift = [ M(:,2:size(M,2)) M(:,size(M,2)) ];
  75. function shift = shiftR(M)
  76.   shift = [ M(:,1) M(:,1:size(M,2)-1) ];
  77. function shift = shiftU(M)
  78.   shift = shiftL(M‘)’;
  79. function S = sussman_sign(D)
  80.   S = D ./ sqrt(D.^2 + 1);
  81. %计算有向距离场函数
  82. function phi = mask2phi(init_a)
  83.   phi=bwdist(init_a)-bwdist(1-init_a)+im2double(init_a);
  84.   phi = –double(phi);
  85. %水平集提取函数,也就是把隐式函数转换成显式函数,所得简单一点,就是提取值为0的等高线
  86. function showCurveAndPhi(I, phi, i)
  87.   imshow(I,‘initialmagnification’,200,‘displayrange’,[ ]);
  88.   hold on;  contour(phi, [0 0], ‘y’,‘LineWidth’,2);
  89.   hold off; title([num2str(i) ‘ Iterations’]); drawnow;
  90. %文献公式2
  91. function f = Dirac2(x, sigma)
  92.     f = (sigma/pi)./(sigma^2+x.^2);
  93. %计算文献公式1
  94. function f = Heaviside2(x, epsilon)
  95.      f = 0.5*(1+2/pi*atan(x./epsilon));
  96. %散度计算
  97. function k = curvature_central(u)
  98. [ux,uy] = gradient(u);
  99. normDu = sqrt(ux.^2+uy.^2+1e-10);
  100. Nx = ux./normDu;
  101. Ny = uy./normDu;
  102. [nxx,junk] = gradient(Nx);
  103. [junk,nyy] = gradient(Ny);
  104. k = nxx+nyy;

 

个人观点:我觉得传统的snake算 法,只能用烂来解释,基本上以遇到一点噪声,就不行了,分割精度真不是一般的差。而水平集的分割只能用高大上来形容,可以进行自动分裂合并等,就是速度有 点慢,因为每次都要对水平集函数进行更新,更新一次就相当于遍历一张图片,因此速度可想而知,当然还有很多加速版的水平集方法,有待测试学习。

转载注明来源:CV视觉网 » 图像处理(九)图像分割(3)泛函能量LevelSet、snake分割

分享到:更多 ()
扫描二维码,给作者 打赏
pay_weixinpay_weixin

请选择你看完该文章的感受:

0不错 0超赞 0无聊 0扯淡 0不解 0路过

评论 4

评论前必须登录!