首先,我见了一个图像处理讨论群,感兴趣的同学们可以加一下,图像处理讨论社: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;