snake算法

能量最小化形式分割的形变模型图像 下面是利用能量最小化形式的形变模型分割 U 形图像的 MATLAB 程序。

[I,map]=rawread('U64.pgm'); %读入 U 形图像 

figure(1); subplot(131); imdisp(I); title('original image');

f=1-I/255; %计算边缘图，使非边缘点灰度值小，边缘点灰度值大 

subplot(132);imdisp(f); title('edge map');

 f0=gaussianBlur(f,1); %与高斯函数做卷积运算，标准差为 1

subplot(133); imdisp(f0); title('gaussianblurred edge map'); %去噪后的边缘图 [px,py]=gradient(f0); %计算高斯外力(梯度)

 figure(2); subplot(121); imdisp(-f); title('gaussian force');subplot(122); 

quiver(px,py); %用 quiver 函数显示箭头图 axis('square', 'equal', 'off', 'ij'); 

figure(2); subplot(121); cla;colormap(gray(64)); image(((1-f)+1)*40); axis('square', 'equal', 'off');

t=0:0.5:6.28; x=32+22*cos(t); y=32+22*sin(t); %给一个初始化圆作为初始曲线

 snakedisp(x,y,'b'); %显示活动轮廓曲线 

[x,y]=snakeinterp(x,y,2,0.5); %进行插值操作，2 为最大间隔，0.5 为最小间隔 

snakedisp(x,y,'r'); %显示活动轮廓曲线

 for i=1:ITER %活动轮廓曲线发生形变的迭代过程

 [x,y]=snakedeform(x,y,1,0.2,10,px,py,5);%调用活动轮廓曲线发生形变子函数 

[x,y]=snakeinterp(x,y,2,0.5); snakedisp(x,y,'r'); title(['Deformation in progress,iter=' num2str(i*5)]);pause(0.5); 

end

figure(3); colormap(gray(64)); image(((1-f)+1)*40);

axis equal snakedisp(x,y,'r'); title(['Final result, iter=' num2str(40*5)]); 

%下面是活动轮廓曲线发生形变的子函数

 function[x,y]=snakedeform(x,y,alpha,beta,gamma,fx,fy,ITER) % alpha 为弹性系数?，beta 为刚性系数?，gamma 为权重因子?，fx 和 fy 为外力场 N=length(x); alpha=alpha*ones(1,N);beta=beta*ones(1,N); % 以下产生 5 对角矢量 

alpham1=[alpha(2:N) alpha(1)]; 

alphap1=[alpha(N) alpha(1:N-1)];

betam1=[beta(2:N) beta(1)]; ]

betap1=[beta(N) beta(1:N-1)]; 

a=betam1; b=-alpha-2*beta-2*betam1; c=alpha+alphap1+betam1+4*beta+betap1; d=-alphap1-2*beta-2*betap1; e=betap1; % 以下产生参数矩阵 

for count=1:ITER  

vfx=interp2(fx,x,y,'*linear'); vfy=interp2(fy,x,y,'*linear');

x=invAI*(gamma*x+vfx); y=invAI*(gamma*y+vfy); 

end

%下面是活动轮廓曲线显示的子函数 functionsnakedisp(x,y,style) 

hold on 

x=x(:); y=y(:); %转为列数据 

if nargin == 3 plot([x;x(1,1)],[y;y(1,1)],style,'MarkerSize',20);

 hold off 

else

disp('snakedisp.m: The input parameter is not correct!'); 

end 

%下面是实现插值操作子函数

function [xi,yi]=snakeinterp(x,y,dmax,dmin) % dmax 为活动轮廓曲线上两个点之间的最大距离，dmin 为最小距离

 x=x(:); y=y(:); N=length(x);d=abs(x([2:N 1])- x(:))+abs(y([2:N 1])- y(:)); %计算活动轮廓曲线上两个点之间的距离 

IDX=(d<dmin); %若距离小于 dmin，则移去一个点

 idx=find(IDX==0); x=x(idx); y=y(idx); N=length(x); d=abs(x([2:N1])- x(:))+abs(y([2:N 1])- y(:)); IDX=(d>dmax); %若距离大于 dmax，则在两个点之间插入一个新的点 z=snakeindex(IDX); %用 snakeindex 函数为活动轮廓曲线上的点编号

 p=1:N+1;xi=interp1(p,[x;x(1)],z'); yi=interp1(p,[y;y(1)],z'); N=length(xi);d=abs(xi([2:N 1])- xi(:))+abs(yi([2:N 1])- yi(:));

while (max(d)>dmax) %在活动轮廓曲线上增加一个新的点

IDX=(d>dmax); z=snakeindex(IDX); p=1:N+1; xi=interp1(p,[xi;xi(1)],z');yi=interp1(p,[yi;yi(1)],z'); N=length(xi); d=abs(xi([2:N 1])-xi(:))+abs(yi([2:N 1])- yi(:)); 

end