图像处理（十一）图像分割(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算法，只能用烂来解释，基本上以遇到一点噪声，就不行了，分割精度真不是一般的差。而水平集的分割只能用高大上来形容，可以进行自动分裂合并等，就是速度有点慢，因为每次都要对水平集函数进行更新，更新一次就相当于遍历一张图片，因此速度可想而知，当然还有很多加速版的水平集方法，有待测试学习。本文地址：/article/7649369.html
作者：hjimce 联系qq：1393852684 更多资源请关注我的博客：http://blog.csdn.net/hjimce
原创文章，版权所有，转载请保留本行信息