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

  • 1
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值