LevalSet水平集分割算法 matlab程序源码

首先，我见了一个图像处理讨论群，感兴趣的同学们可以加一下，图像处理讨论社：397489395

基于水平集的图像分割算法算是一种进化的snake算法，是一种隐式函数的方法，目前代码以及原理还处于深究阶段，有些问题没有搞懂，代码摘自 https://blog.csdn.net/hjimce/article/details/45586727；欢迎大家一同探讨图像处理问题

function seg = local_AC_MS(Img, mask_init, rad, alpha, num_it, epsilon )
  %Img为灰度图像；mask_init为初始曲线的一个mask,曲线上的点坐标在mask中的值为1；
  %rad为卷积半径；alpha和epsilon为公式中参数， num_it为迭代次数
  
  phi0 = -double(bwdist(mask_init)-bwdist(1-mask_init)+im2double(mask_init));  %计算有限距离场，将显示曲线转化为水平集函数
  phi = phi0;
  
B0 = ones(2*rad+1, 2*rad+1); % 卷积核mask
KI = conv2(Img, B0, 'same');  %对Img图像进行卷积
KONE = conv2(ones(size(Img), B0, 'same'));

%计算使曲线进行演化的公式
for ii = 1:num_it
    mask = 0.5*(1+2/pi*atan(phi./epsilon));% 计算封闭曲线的内外部
    
    I = Img.*mask;%区分曲线的内外部分
    temp1 = conv2(mask,B0,'same');
    temp2 = conv2(I, B0, 'same');
    c1 = temp2./(temp1);
    c2 = (KI-temp2)./(KONE-temp1);
    
    A1 = temp1;
    A2 = conv2(1-mask,B0,'same');
    D = (A1.*A2+eps);  %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');
    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);
    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;
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 = zreos(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;
dt .* sussman_sign(D) .* dD;

function shift = shiftD(M)
shiftR(M')';

function shift = shiftL(M)
[M(:,2:size(M,2)) M(:,size(M,2))];

function shift = shiftR(M)
[M(:,1) M(:,1:size(M,2)-1)];

function shift = shiftU(M)
shiftL(M');

function S = sussman_sign(D)
S = D ./ sqrt(D.^2 + 1);

%水平集提取函数，也就是把隐式函数转换为显示函数，所得简单一点，就是提取值为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;

%散度计算
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;