能量最小化形式分割的形变模型图像 下面是利用能量最小化形式的形变模型分割 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