天津大学研究生数字图像处理作业-Snake模型..

Snake模型简介及其编程实现

Snake模型也称为主动轮廓线模型，最初由Kass等人在1987年第一届计算机国际视觉会议上提出，一经提出就成为计算机视觉领域研究的热点。Snake的基本思想是通过人的识别能力，在图像中目标边界附近确定初始轮廓线，然后对曲线进行能量最小化变形，使其锁定在待分割目标的边界上。Snake模型之所以能得到如此重视，是因为它将图像目标的先验知识（如大小、位置、形状等）与图像特征（灰度、梯度、纹理等）结合起来，克服了传统图像分割方法将二者分离的缺陷。近年来，许多文章从传统Snake模型的能量函数构造和求解算法方面进行改进，在其基础上衍生出了许多新的Snake模型。

1、Snake模型的基本原理

其基本思想是依据图像信息进行曲线（曲面）演化，使其最终找到目标物体的边界。这种方法将分割问题转化为最优化问题，利用闭合曲线（或曲面）形变的特定规律，定义度量闭合曲线（曲面）形变的能量函数，通过最小化能量函数使曲线（曲面）逐渐逼近图像中目标物体的边缘。先提供待分割图像的一个初始轮廓的位置，并对其定义个能量函数，是轮廓沿能量降低的方向靠近。当能量函数达到最小的时候，提供的初始轮廓收敛到图形中目标的真实轮廓。

Snake能量函数是有内部能量函数和外部能量函数组成，内部能量控制轮廓的平滑性和连续性，外部能量由图像能量和约束能量组成，控制轮廓向着实际轮廓收敛，其中约束能量可根据具体的对象形态定义，使得snake具有很大的灵活性。

Snake模型发展10多年来，许多学者对于经典的snake模型做了改进，提出各种改进的snake模型，其中梯度矢量流（Gradient Vector Flow, GVF）模型扩大了经典snake的外力作用范围，加强了对目标凹轮廓边缘的吸引力，提高了传统的snake模型。

2. 基本的Snake模型

数学上，将活动轮廓表示成一条参数曲线V(s,t)=(x(s,t),y(s,t))，其中，V是曲线点的二维坐标，t是时间参数，s是弧长参数。轮廓的总弧长归一化到1。改曲线的能量可以用能量泛函表示为E(V)=E int(V)+E ext(V)，E int(V)是内部能量泛函，E ext(V)是外部能量泛函。曲线V在图像的空间域运行使得E(V)最小。

其中内部能量泛函定义为：

两个物理参量表明曲线的物理特征：

是张紧系数，值越大，轮廓曲线收缩越快；是强度系数，值越大，轮廓曲线变得僵硬而不易弯曲。

外部能量泛函定义为：

其中E ext (V)是定义在图像平面上的标量势函数，设I(x,y)是一灰度图像，一个可能使活动轮廓想边界运动的外部能量函数可以设计

为

,其中，c 是一个正的权因子，它控制势的

大小，

是梯度算子，是标准差为的二维高斯函数，*是二维图像卷积算子。

表明图像进过一个高斯平滑滤波，用它代替I 可降低计算梯度的噪

声。

根据变分原理使能量泛函E(V)最小化的曲线V(S)满足下述欧拉方程： 该方程可解释为力平衡方程，它表示当轮廓达到平衡点时内力和外力的平衡。其中，前两项参数表示内部的拉伸弹性力和弯曲刚性力，第三项参数表示曲线所受到的图像信息外力。上式表明曲线能量最小化过程就是在图像信息外力和曲线本身的内力作用下运动，达到平衡状态的过程。

3. 改进的Snake 模型

传统的snake 模型存在不足之处：首先，初始轮廓线必须接近真实的图像边缘。其次，Snake 模型有可能收敛到局部极值点，无法逼近物体内部的“凹形”轮廓。针对这些缺陷，研究者主要从能量函数构造和求解方面进行改进。

（1）. 能量函数构造方面的改进。能量函数构造方面的改进具有代表性的主要有Cohen 等人提出的气球Snake 模型和Xu 提出的GVF Snake 模型。气球Snake 模型在Snake 模型中增加了气球力，大小为常数，方向沿轮廓点的法线外方向。在气球力的作用下，轮廓线作为一个整体膨胀或收缩，当轮廓线进入图像能量场的作用范围后，被吸引向感兴趣区域的边界。Cohen 的气球力模型力场吸引范围大，克服了初始化缺陷，但为了使克服噪声和不越过弱边界两者达到平衡，在怎样选取合适的气球力权值方面仍很难把握。Cohen 等人还提出一种距离势能0ss ssss image v v E α-β-?

=

模型，虽然改进了初始化问题，但Snake曲线无法收敛到凹陷区域。Xu采用梯度向量场作为模型的外部能量场，产生了GVF Snake模型该模型扩大了Snake 的捕获区，初始化轮廓可以不必充分靠近真实边缘，并能将Snake模型轮廓曲线拖向物体的深度凹陷区。但是该算法计算量太大，实时性较差。

（2）. 数值计算方法的改进。Kass采用变分法解能量函数方程，一般只能得到局部最优解，还可能出现数值解不稳定的现象，并且要求外力必须是可微的，使算法的应用受到限制。Amini采用动态规划算法来求解局部最优解，改方法数值稳定，还可以增加硬约束，放宽了对外力可微的限制，但运算速度慢，时间复杂度高。Williams在此基础上应用贪婪算法，在保持解的稳定性情况下，将时间复杂度明显降低。贪婪算法实现简单，速度快，容易加入新约束，最终结点分布比较均匀。贪婪算法虽然有上述优点，但在应用中仍然存在一些问题：没有明确硬强制力对Snake模型的作用方法，也没有给出具体的硬强制力模型；代表轮廓线的控制点点数固定不变，没有给出初始轮廓点的选取方法。

目前，许多有关Snake模型的文章主要是在上述几种模型的基础上进行改进，采用较多的是Xu提出的梯度向量流Snake模型和Williams的贪婪算法。4. Snake的进化模型

(1). McInerney 提出一种拓扑自适应snake模型（Topology Adaptive Snake,T-Snake）

该算法基于仿射细胞图像分解（Affine Cell Image Decomposition，ACID）先在待分割图像上加以三角网格，然后在图像区域的适当位置做一条初始曲线，最后取曲线与网格的交点作为snake的初始离散点，其第i个snake的离散点的坐标为其中，相邻两点，之间由一条弹性样条连接而成

由于T-Snake模型可借助三角形网格和网格点的特征函数来确定边界三角形，可促使snake模型演化过程中的分裂和合并，从而保证了其具有能够处理拓扑结果复杂图像的能力，因此能够很好的满足医学图像拓扑结果复杂的特点。此算法用于脑部MR切片有良好的性能。

(2). 双T-Snake模型

双T-Snake模型（Dual-T-Snakes）是在T-Snake模型的基础上产生的，其主要思想是采用内外两个初始轮廓，其中一个轮廓从目标外向内收缩和分裂，另一个轮廓从目标内部向外膨胀，两个初始轮廓可以离目标边界较远，迭代的过程中对能量较大的轮廓增加驱动力，使其靠近与之相对应的轮廓，直到连个轮廓收敛到同一个为止

(3). Loop Snake 模型

Loop Snake模型是一种加强了拓扑控制的T-Snake模型，这种方法的关键集中在曲线的每一步进化中都要形成循环，其基本思想是，确保图像轮廓曲线精

确地线性地映射到适当的分类中，然后在额外的记过loop-Tree的帮助下，尽可能少的时间内运用已经被snake探究的循环来决定是否进行区域划分，这种模型的实质是对T-Snake模型的一种改进。由于加强了拓扑控制，使得Loop Snake 模型既可以忽略背景中强噪声又可以在演化过程中进行多次分裂。

(4). 连续snake模型

在Snake模型中，轮廓曲线由一条给定容许误差范围的光滑曲线组成，相对于离散snake来说，连续snake模型所需要的控制点少，比离散的更具优越性。

(5). B-Snake模型

B-Snake模型是通过B样条曲线来定义的，其轮廓曲线由各曲线段光滑相连而成，每一个曲线段都是由一个给定次数多项式表示，这种多项式是B样条曲度函数的一种线性组合，并以控制点为系数。在有些B-Snake模型中并没有明确应用内部能量，这是因为B样条本身就含有内部能量，snake轮廓曲线只受外力影响着图像边缘移动。可用于对图像切片分割区域的描述与跟踪而用于器官的三维重建。

5. 基本Snake模型的程序实现

使用matlab实现基本的snake模型，主程序如下（主程序中调用的部分函数将在附录中给出）：

cd ..; s = cd; s = [s, '/snake']; path(s, path); cd examples;

help tradition_ex;

% ==== Example 1: U-shape object ====

% 读取64x64 U形图像

[I,map] = rawread('../images/U64.pgm');

%边缘检测

disp(' 边缘检测...');

f = 1 - I/255; %图像归一化

f0 = gaussianBlur(f,1); %高斯函数模糊处理

% note: snake potential is the negative of edge map

disp(' 计算外部力...');

[px,py] = gradient(f0);%其中px为其水平方向上的梯度，py为其垂直方向上的梯度

% 显示结果

figure(1);

subplot(121); imdisp(-f); title('snake potential');

subplot(122); quiver(px,py); %在x-y平面上绘制向量（px,py）axis('image', 'off', 'ij'); % 确定图像显示坐标

title('traditional force');

% snake变换

disp(' ');

disp(' 请按任意键开始变换：');

pause;

figure(1); subplot(121); %将图一分为二，并在左侧显示一张。colormap(gray(64));%得出一个64色的灰度图

image(((1-f)+1)*40); %产生正方形坐标系

axis('square', 'off');

disp('');

disp('... 捕获范围正在缩小...')

t = 0:0.05:6.28;

x = 32 + 17*cos(t);

y = 32 + 17*sin(t);

[x,y] = snakeinterp(x,y,3,1); %插入或移动节点处理snakedisp(x,y,'r')

pause(1);

disp('... 按-C 终止程序！');

for i=1:100,

[x,y] = snakedeform(x,y,0.05,0,1,4,px,py,5);

[x,y] = snakeinterp(x,y,3,1);

snakedisp(x,y,'r')

title(['变换中..., 迭代次数= ' num2str(i*5)])

pause(0.5);

end

disp(' ');

disp(' 请按任意键显示最终结果：');

pause;

subplot(121);

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

axis('square', 'off');

snakedisp(x,y,'r')

title(['最终结果：迭代次数= ' num2str(i*5)]); 程序运行结果如下：

图1. 原始图像及图像的传统力场

图2. Snake模型变换中

图3. Snake模型变换的最终结果

由实验结果可以看出，传统高斯力Snake模型在迭代求解过程中，很好的使初始轮廓线向U形图的边缘靠近，最终得到U形图的边缘轮廓；同时，从图中可以看出，可变轮廓线无法捕捉U形的底部，轮廓没有收缩到U形图的凹陷部分。

附录：

1.rawread.m

function [X,map] = rawread(filename,n,m)

% RAWREAD 读入可移植的点阵图文件, 或者原始文件.

%RAWREAD('imagefile.raw', xsize, ysize) 读入"raw"格式的图像文件

%RAWREAD('imagefile.pgm') 读入一"pgm" 格式的图像文件

%[X,map] = RAWREAD('imagefile.raw') 同时返回一个灰度图和一个彩色图%[X,map] = rawread('imagefile.raw',sx,sy);

%or [X,map] = rawread('imagefile.pgm');

% image(X)

% colormap(map) 显示出正确的色彩。

dot = max(find(filename == '.'));

suffix = filename(dot+1:dot+3);

if strcmp(suffix,'pgm') | strcmp(suffix,'raw')%strcmp比较字符

disp(sprintf('\nopens %s file\n',filename));

fp = fopen(filename,'rb','b'); % "Big-endian" byte order.

if (fp<0)

error(['Cannot open ' filename '.']);

end

if strcmp(suffix,'pgm')

% Read and crack the header(打开标题）

head = fread(fp,2,'uchar'); % pgm magic number : P5

if ~strcmp(head,'P5'),

fprintf(1,'\n Magic Number : %s\n',head);

else

fprintf(1,'\n Bad Magic Number : %s\n',head);

error('cannot continue this way, good bye cruel world');

end

c = fread(fp,1,'uchar'); %reads the carriage return（回车符）separating P5 from the creator

precreator = fread(fp,1,'uchar'); % look for a '#' character preceeding a creator signature

if precreator == '#',

c = setstr(20); % any character except carriage return

cr = setstr(10); % defines a carriage return

while c ~= cr,

c = fread(fp,1,'uchar');

creator = [creator,c];

end;

fprintf(1,'\n creator : %s\n',creator);

else

fprintf('\n No creator signature\n');

fseek(fp,-1,'cof'); % return one char before

end;

end

if nargin <2,

if strcmp(suffix,'raw')

% assume image size is 256x256

disp('RAW file without size : assume image size is 256x256');

n = 256;

m = 256;

else % for PGM files

%读入大小和深度

disp(' reads sizes');

n = fscanf(fp,'%d',1);

tn = num2str(n);

disp([' xsize = ' tn]);

m = fscanf(fp,'%d',1);

tm = num2str(m);

disp([' ysize = ' tm]);

p = fscanf(fp,'%d',1);

tp = num2str(p);

disp([' depth = ' tp]);

c = fread(fp,1,'uchar'); %reads the last carriage return

end;

end

%生成一个灰度调色板，区间为[0,1].

disp(' create gray palette');

for i=1:256,

map(i,[1:3])=[i/256,i/256,i/256];

end;

% 读取图像

disp(' Reads image data ...');

[X,l] = fread(fp,[n,m],'uchar');

if l ~= m*n, l, error('HSI image file is wrong length'), end

X = X'+1;

fclose(fp);

disp('end');

else

error('Image file name must end in ''raw'' or ''pgm''.')

end

2.gaussianBlur.m

function GI = gaussianBlur(I,s)

% GAUSSIANBLUR 用高斯核函数模糊处理图像

% GI = gaussianBlur(I,s)中I是图像

% s 是高斯核函数的标准误差值，GI 是高斯模糊处理后的图像。

M = gaussianMask(1,s);%高斯掩膜处理。

M = M/sum(sum(M)); % 规格化高斯掩膜，使和为1.

GI = xconv2(I,M);

3.snakeinterp.m

function [xi,yi] = snakeinterp(x,y,dmax,dmin)

% [xi,yi] = snakeinterp(x,y,dmax,dmin)

% dmax: snake点之间的最大距离。

% dmin: snake点之间的最小距离。

% 如果d(i,i+1)>dmax, 则要在i与i+1之间添加一个点。

% 如果d(i,i+1)

x = x(:); y = y(:); %分别将x,y以列向量的形式输出

N = length(x);

d = abs(x([2:N 1])- x(:)) + abs(y([2:N 1])- y(:));

% remove the points which distance to neighbor points is shorter than dmin IDX = (d

idx = find(IDX==0);%实际上idx=1:126（列向量）

x = x(idx);

y = y(idx);

N = length(x);

d = abs(x([2:N 1])- x(:)) + abs(y([2:N 1])- y(:));

IDX = (d>dmax);

z = snakeindex(IDX);%本例中z=1:N,输出z=1:N

p = 1:N+1;

xi = interp1(p,[x;x(1)],z');%interp1(库函数）；这里[x;x(1)]的表述是对的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

4.snakedisp.m

function snakedisp(x,y,style)

hold on

% 转换为列向量

x = x(:); y = y(:);

if nargin == 3

plot([x;x(1,1)],[y;y(1,1)],style);

hold off

else

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

end

5.snakedeform.m

function [x,y] = snakedeform(x,y,alpha,beta,gamma,kappa,fx,fy,ITER) generates the parameters for snake

% alpha=0.05 beta=0.1 gamma=0.5 ITER=5

N = length(x);

alpha = alpha* ones(1,N);

beta = beta*ones(1,N);

%产生对角向量

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;

% 生成参数矩阵

A = diag(a(1:N-2),-2) + diag(a(N-1:N),N-2);

A = A + diag(b(1:N-1),-1) + diag(b(N), N-1);

A = A + diag(c);

A = A + diag(d(1:N-1),1) + diag(d(N),-(N-1));

A = A + diag(e(1:N-2),2) + diag(e(N-1:N),-(N-2));

invAI = inv(A + gamma * diag(ones(1,N)));

for count = 1:ITER,

vfx = interp2(fx,x,y,'*linear'); %interpolation（插入）

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

% snake变换

x = invAI * (gamma* x + kappa*vfx);

y = invAI * (gamma* y + kappa*vfy);

end