2014-06-06 16:21:13 gujinjin2008 阅读数 17602
% =========================================================================
%                   Snakes:Active Contour Models
% =========================================================================
% By gujinjin 2012/12/10-12  Sunny
% 基于KASS等的论文思想
% 参考文献:
% [1] KASS etc.Snakes:Active Contour Models
% [2] CSDN 博客 - Author:乐不思蜀Tone
% [3] Ritwik Kumar(Harvard University),D.Kroon(Twente University)的程序
% [4] 《数学建模与数学实验》
% ------
clc;clf;clear all;

% =========================================================================
%                      获取手动取点坐标
% =========================================================================
% 读取显示图像
%I = imread('Coronary_MPR.jpg');
I =  imread('test1.png');
% 转化为双精度型
%I = im2double(I); 
% 若为彩色,转化为灰度
%if(size(I,3)==3), I=rgb2gray(I); end
figure(1),imshow(I);
%---------------------------
%        高斯滤波
%---------------------------
sigma=1;
% 创建特定形式的二维高斯滤波器H
H = fspecial('gaussian',ceil(3*sigma), sigma);
% 对图像进行高斯滤波,返回和I等大小矩阵
Igs = filter2(H,I,'same');
%---------------------------
%      获取Snake的点坐标
%---------------------------
figure(2),imshow(Igs);
x=[];y=[];c=1;N=100; %定义取点个数c,上限N
% 获取User手动取点的坐标
% [x,y]=getpts
while c<N
    [xi,yi,button]=ginput(1);
    % 获取坐标向量
    x=[x xi];
    y=[y yi];
    hold on
    % text(xi,yi,'o','FontSize',10,'Color','red');
    plot(xi,yi,'ro');
    % 若为右击,则停止循环
    if(button==3), break; end
    c=c+1;
end

% 将第一个点复制到矩阵最后,构成Snake环
xy = [x;y];
c=c+1;
xy(:,c)=xy(:,1);
% 样条曲线差值
t=1:c;
ts = 1:0.1:c;
xys = spline(t,xy,ts);
xs = xys(1,:);
ys = xys(2,:);
% 样条差值效果
hold on
temp=plot(x(1),y(1),'ro',xs,ys,'b.');
legend(temp,'原点','插值点');

% =========================================================================
%                     Snakes算法实现部分
% =========================================================================
NIter =1000; % 迭代次数
alpha=0.2; beta=0.2; gamma = 1; kappa = 0.1;
wl = 0; we=0.4; wt=0;
[row col] = size(Igs);

% 图像力-线函数
Eline = Igs;
% 图像力-边函数
[gx,gy]=gradient(Igs);
Eedge = -1*sqrt((gx.*gx+gy.*gy));
% 图像力-终点函数
% 卷积是为了求解偏导数,而离散点的偏导即差分求解
m1 = [-1 1]; 
m2 = [-1;1];
m3 = [1 -2 1]; 
m4 = [1;-2;1];
m5 = [1 -1;-1 1];
cx = conv2(Igs,m1,'same');
cy = conv2(Igs,m2,'same');
cxx = conv2(Igs,m3,'same');
cyy = conv2(Igs,m4,'same');
cxy = conv2(Igs,m5,'same');

for i = 1:row
    for j= 1:col
        Eterm(i,j) = (cyy(i,j)*cx(i,j)*cx(i,j) -2 *cxy(i,j)*cx(i,j)*cy(i,j) + cxx(i,j)*cy(i,j)*cy(i,j))/((1+cx(i,j)*cx(i,j) + cy(i,j)*cy(i,j))^1.5);
    end
end

%figure(3),imshow(Eterm);
%figure(4),imshow(abs(Eedge));
% 外部力 Eext = Eimage + Econ
Eext = wl*Eline + we*Eedge + wt*Eterm;
% 计算梯度
[fx,fy]=gradient(Eext);

xs=xs';
ys=ys';
[m n] = size(xs);
[mm nn] = size(fx);

% 计算五对角状矩阵
% 附录: 公式(14) b(i)表示vi系数(i=i-2 到 i+2)
b(1)=beta;
b(2)=-(alpha + 4*beta);
b(3)=(2*alpha + 6 *beta);
b(4)=b(2);
b(5)=b(1);

A=b(1)*circshift(eye(m),2);
A=A+b(2)*circshift(eye(m),1);
A=A+b(3)*circshift(eye(m),0);
A=A+b(4)*circshift(eye(m),-1);
A=A+b(5)*circshift(eye(m),-2);

% 计算矩阵的逆
[L U] = lu(A + gamma.* eye(m));
Ainv = inv(U) * inv(L); 

% =========================================================================
%                      画图部分
% =========================================================================
%text(10,10,'+','FontName','Time','FontSize',20,'Color','red');
% 迭代计算
figure(3)
for i=1:NIter;
    ssx = gamma*xs - kappa*interp2(fx,xs,ys);
    ssy = gamma*ys - kappa*interp2(fy,xs,ys);
 
    % 计算snake的新位置
    xs = Ainv * ssx;
    ys = Ainv * ssy;
    
    % 显示snake的新位置
    imshow(I); 
    hold on;
    plot([xs; xs(1)], [ys; ys(1)], 'r-');
    hold off;
    pause(0.001)    
end



2010-10-15 12:57:00 BeanYoung 阅读数 2935

最近在实现GVF snake,其中要对读取的图像求解GVF力场。图像以灰度模式读进来之后(IPL_DEPTH_8U),要进行归一化处理,转化成深度为IPL_DEPTH_32F的图像,但是现实不正常。google到了这篇文章 OpenCV对不同图像深度的处理 ,文章中对opencv里面的不同图像深度进行了分析和测试。如下:

 

      测试double型:0.0--1.0之间                           IPL_DEPTH_64F
      测试float型:0.0--1.0之间                              IPL_DEPTH_32F
      测试long型:0--65535之间                             IPL_DEPTH_32S       
      测试short int型:-32768--32767之间                  IPL_DEPTH_16S       
      测试unsigned short int型:0--65535之间              IPL_DEPTH_16U
      测试char型:-128--127之间                            IPL_DEPTH_8S         
      测试unsigned char型:0--255之间                     IPL_DEPTH_8U

 

在进行不同深度图像间进行转换的时候,如果只是调用cvScale()函数的话,用cvShowImage()来显示图像会不正常。调用过cvScale()之后,要将转换过后的图像进行处理,才能正常的用cvShowImage()函数显示。如下代码显示了从IPL_DEPTH_8U转换到IPL_DEPTH_32F后显示。经归一化处理后便可正常显示。也可用其他方法进行处理,这里使用归一化是因为后面要用到图像的归一化,我这里就顺便贴了归一化的代码。

 

ps: CSDN上下载的那个gvf snake c++的代码有点问题,根本不能收敛到凹区域,增加迭代次数后也不能收敛到凹区域。而且CvSepFilter在2.1版本中也找不到,用cvFilter2D()替换了CvSepFilter,问题仍不能得到解决。且cvSobel()处理edge map时会有溢出。

预计GVF snake最近一个月能够实现。

 

小抱怨一下,这学期课多,导师催的也紧,图像处理也是刚刚接触。能逃的课都逃了,逃课都成了习惯,不知最后期末考试成绩会是个什么样。

 

 

2017-06-12 22:29:15 coming_is_winter 阅读数 9220

图像处理之图像分割(一)之活动轮廓模型:Snake算法简单梳理


  Snake算法,应该也可以翻译成蛇形算法,或者是包含曲折前进的意思。具体函数背景原理介绍参考:zouxy09,http://blog.csdn.net/zouxy09/article/details/8712287,图像分割之(五)活动轮廓模型之Snake模型简介,这里还是就几个自己思考的重点拓展一下。

  

 1. 这段中讲v(s)的表示,s∈[0,1],和s是以傅里叶变换形式描述边界的自变量。自己理解不是很明白,查了点资料,整合一下,说法有可能不严谨。首先:




   2.Snake模型到底是怎么变化的?



  弹性力量和弯曲力量合成内部力量,,保持固有的自身形状,外部力量改变其固有的自身形状,在合适的范围内受某种力而变化,这里举个例子:
 


  上面图中这种特大气泡(或想一下气球)的形状是有什么决定的呢?其中的弹性力量是曲线的斜率使其保持大概的圆形(可以理解为材料的张力,失重条件下水会变成水球),虽然现在上图的圆形变形严重,而弯曲力量是曲线斜率的斜率,这里理解为收缩的力量,可以想象一下气球紧绷的力量(这是气球本身材质决定,与外力无关),如果有漏气现象,紧绷的力量会使气泡或气球逐渐收缩,直至最小。外部力量这里的是空气的压力,这在上图中可以看到气泡的不规则凹凸,可以说是大气的压力的力量,形象的说是空气力量梯度的变化,在Snake中表现为像素灰度值梯度的变化。联系再丰富一点,Snake迭代的过程可以想一下真空包装:



  一开始先验知识人为划定Snake框架,相当于还未抽气时的包装袋,抽气开始在内力(包装袋的形状与材质决定)与外力(大气压力)的作用下,不断迭代,最终的轮廓收缩到鸡腿形状,能量最小达到平衡。

  3.泛函与变分
  函数的函数为泛函,我们求解的能量函数公式最终表现为函数的函数,怎样求公式最小值,一般求导找零点。函数的函数求导为变分。

  参考文献:
  [数字图像处理(第二版)].(美)冈萨雷斯.扫描版。

2019-09-21 21:10:20 webzhuce 阅读数 510

示例演示

        如果在中文搜索的话,一般会找到《数字图像处理-图像分割:Snake主动轮廓模型 Matlab代码及运行结果》。里面有句代码,千万别用,否则出不来效果。(别问我怎么知道的)

% 转化为双精度型
%I = im2double(I);

         当然英文搜索的话,会在Matlab官网上找到《Snakes: Active Contour Models》,提供更好的代码,里面不同图片的算法参数都调好了。核心代码展示。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% interate.m implements the core of the snakes (active contours) algorithm.
% It is called by the snk.m file which is the GUI frontend. If you do not
% want to deal with GUI, look at this file and the included comments which
% explain how active contours work.
%
% To run this code with GUI
%   1. Type guide on the matlab prompt.
%   2. Click on "Go to Existing GUI"
%   3. Select the snk.fig file in the same directory as this file
%   4. Click the green arrow at the top to launch the GUI
%
%   Once the GUI has been launched, you can use snakes by
%   1. Click on "New Image" and load an input image. Samples image are
%   provided.
%   2. Set the smoothing parameter "sigma" or leave it at its default value
%   and click "Filter". This will smooth the image.
%   3. As soon as you click "Filter", cross hairs would appear and using
%   them and left click of you mouse you can pick initial contour location
%   on the image. A red circle would appead everywhere you click and in
%   most cases you should click all the way around the object you want to
%   segement. The last point must be picked using a right-click in order to
%   stop matlab for asking for more points.
%   4. Set the various snake parameters (relative weights of energy terms
%   in the snake objective function) or leave them with their default value
%   and click "Iterate" button. The snake would appear and move as it
%   converges to its low energy state.
%
% Copyright (c) Ritwik Kumar, Harvard University 2010
%               www.seas.harvard.edu/~rkkumar
%
% This code implements 揝nakes: Active Contour Models�by Kass, Witkin and
% Terzopolous incorporating Eline, Eedge and Eterm energy factors. See the
% included report and the paper to learn more.
%
% If you find this useful, also look at Radon-Like Features based
% segmentation in  the following paper:
% Ritwik Kumar, Amelio V. Reina & Hanspeter Pfister, 揜adon-Like Features 
% and their Application to Connectomics� IEEE Computer Society Workshop %
% on Mathematical Methods in Biomedical Image Analysis (MMBIA) 2010
% http://seas.harvard.edu/~rkkumar
% Its code is also available on MATLAB Central
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [smth] = interate(image, xs, ys, alpha, beta, gamma, kappa, wl, we, wt, iterations);
% image: This is the image data
% xs, ys: The initial snake coordinates
% alpha: Controls tension
% beta: Controls rigidity
% gamma: Step size
% kappa: Controls enegry term
% wl, we, wt: Weights for line, edge and terminal enegy components
% iterations: No. of iteration for which snake is to be moved


%parameters
N = iterations; 
smth = image;

% Calculating size of image
[row col] = size(image);


%Computing external forces

eline = smth; %eline is simply the image intensities

[grady,gradx] = gradient(smth);
eedge = -1 * sqrt ((gradx .* gradx + grady .* grady)); %eedge is measured by gradient in the image

%masks for taking various derivatives
m1 = [-1 1];
m2 = [-1;1];
m3 = [1 -2 1];
m4 = [1;-2;1];
m5 = [1 -1;-1 1];

cx = conv2(smth,m1,'same');
cy = conv2(smth,m2,'same');
cxx = conv2(smth,m3,'same');
cyy = conv2(smth,m4,'same');
cxy = conv2(smth,m5,'same');

for i = 1:row
    for j= 1:col
        % eterm as deined in Kass et al Snakes paper
        eterm(i,j) = (cyy(i,j)*cx(i,j)*cx(i,j) -2 *cxy(i,j)*cx(i,j)*cy(i,j) + cxx(i,j)*cy(i,j)*cy(i,j))/((1+cx(i,j)*cx(i,j) + cy(i,j)*cy(i,j))^1.5);
    end
end

% imview(eterm);
% imview(abs(eedge));
eext = (wl*eline + we*eedge -wt * eterm); %eext as a weighted sum of eline, eedge and eterm

[fx, fy] = gradient(eext); %computing the gradient


%initializing the snake
xs=xs';
ys=ys';
[m n] = size(xs);
[mm nn] = size(fx);
    
%populating the penta diagonal matrix
A = zeros(m,m);
b = [(2*alpha + 6 *beta) -(alpha + 4*beta) beta];
brow = zeros(1,m);
brow(1,1:3) = brow(1,1:3) + b;
brow(1,m-1:m) = brow(1,m-1:m) + [beta -(alpha + 4*beta)]; % populating a template row
for i=1:m
    A(i,:) = brow;
    brow = circshift(brow',1)'; % Template row being rotated to egenrate different rows in pentadiagonal matrix
end

[L U] = lu(A + gamma .* eye(m,m));
Ainv = inv(U) * inv(L); % Computing Ainv using LU factorization
 
%moving the snake in each iteration
for i=1:N;
    
    ssx = gamma*xs - kappa*interp2(fx,xs,ys);
    ssy = gamma*ys - kappa*interp2(fy,xs,ys);
    
    %calculating the new position of snake
    xs = Ainv * ssx;
    ys = Ainv * ssy;
    
    
    %Displaying the snake in its new position
    imshow(image,[]); 
    hold on;
    
    plot([xs; xs(1)], [ys; ys(1)], 'r-');
    hold off;
    pause(0.001)    
end;

不同图片运行结果
在这里插入图片描述
在这里插入图片描述

在这里插入图片描述

2016-09-13 10:16:47 sinat_24206709 阅读数 12657

《Matlab图像处理》part1 Snakes:Active Contour Models 主动轮廓模型

参考博客:

数字图像处理-图像分割:Snake主动轮廓模型 Matlab代码及运行结果

图像分割之(五)活动轮廓模型之Snake模型简介


简介

在“图像分割之(一)概述”中咱们简单了解了目前主流的图像分割方法。下面咱们主要学习下基于能量泛函的分割方法。这里学习下Snake模型简单的知识,Level Set(水平集)模型会在后面的博文中说到。

基于能量泛函的分割方法:

该类方法主要指的是活动轮廓模型(active contour model)以及在其基础上发展出来的算法,其基本思想是使用连续曲线来表达目标边缘,并定义一个能量泛函使得其自变量包括边缘曲线,因此分割过程就转变为求解能量泛函的最小值的过程,一般可通过求解函数对应的欧拉(Euler.Lagrange)方程来实现,能量达到最小时的曲线位置就是目标的轮廓所在。

主动轮廓线模型是一个自顶向下定位图像特征的机制,用户或其他自动处理过程通过事先在感兴趣目标附近放置一个初始轮廓线,在内部能量(内力)和外部能量(外力)的作用下变形外部能量吸引活动轮廓朝物体边缘运动,而内部能量保持活动轮廓的光滑性和拓扑性,当能量达到最小时,活动轮廓收敛到所要检测的物体边缘。

一、曲线演化理论


曲线演化理论在水平集中运用到,但我感觉在主动轮廓线模型的分割方法中,这个知识是公用的,所以这里我们简单了解下。

曲线可以简单的分为几种:


曲线存在曲率,曲率有正有负,于是在法向曲率力的推动下,曲线的运动方向之间有所不同:有些部分朝外扩展,而有些部分则朝内运动。这种情形如下图所示。图中蓝色箭头处的曲率为负,而绿色箭头处的曲率为正。

简单曲线在曲率力(也就是曲线的二次导数)的驱动下演化所具有的一种非常特殊的数学性质是:一切简单曲线,无论被扭曲得多么严重,只要还是一种简单曲线,那么在曲率力的推动下最终将退化成一个圆,然后消逝(可以想象下,圆的所有点的曲率力都向着圆心,所以它将慢慢缩小,以致最后消逝)。
描述曲线几何特征的两个重要参数是单位法矢和曲率,单位法矢描述曲线的方向,曲率则表述曲线弯曲的程度。曲线演化理论就是仅利用曲线的单位法矢和曲率等几何参数来研究曲线随时间的变形。曲线的演变过程可以认为是表示曲线在作用力 F 的驱动下,朝法线方向 N 以速度 v 演化。而速度是有正负之分的,所以就有如果速度 v 的符号为负,表示活动轮廓演化过程是朝外部方向的,如为正,则表示朝内部方向演化,活动曲线是单方向演化的,不可能同时往两个方向演化。

所以曲线的演变过程,就是不同力在曲线上的作用过程,力也可以表达为能量。世界万物都趋向于能量最小而存在。因为此时它是最平衡的,消耗最小的(不知理解对不?)。那么在图像分割里面,我们目标是把目标的轮廓找到,那么在目标的轮廓这个地方,整个轮廓的能量是最小的,那么曲线在图像任何一个地方,都可以因为力朝着这个能量最小的轮廓演变,当演变到目标的轮廓的时候,因为能量最小,力平衡了,速度为0了,也就不动了,这时候目标就被我们分割出来了。

那现在关键就在于:1)这个轮廓我们怎么表示;2)这些力怎么构造,构造哪些力才可以让目标轮廓这个地方的能量最小?

这两个问题的描述和解决就衍生出了很多的基于主动轮廓线模型的分割方法。第一个问题的回答,就形成了两大流派:如果这个轮廓是参数表示的,那么就是参数活动轮廓模型(parametric active contour model),典型为snake模型,如果这个轮廓是几何表示的,那么就是几何活动轮廓模型(geometric active contour model),即水平集方法(Level Set),它是把二维的轮廓嵌入到三维的曲面的零水平面来表达的(可以理解为一座山峰的等高线,某个等高线把山峰切了,这个高度山峰的水平形状就出来了,也就是轮廓了),所以低维的演化曲线或曲面,表达为高维函数曲面的零水平集的间接表达形式(这个轮廓的变化,直观上我们就可以调整山峰的形状或者调整登高线的高度来得到)。

那对于第二个问题,是两大流派都遇到的问题,是他们都需要解决的最关键的问题。哪些力才可以达到分割的目标呢?这将在后面聊到。


二、Snakes模型


自1987年Kass提出Snakes模型以来,各种基于主动轮廓线的图像分割理解和识别方法如雨后春笋般蓬勃发展起来。Snakes模型的基本思想很简单,它以构成一定形状的一些控制点为模板(轮廓线),通过模板自身的弹性形变,与图像局部特征相匹配达到调和,即某种能量函数极小化,完成对图像的分割。再通过对模板的进一步分析而实现图像的理解和识别。

简单的来讲,SNAKE模型就是一条可变形的参数曲线及相应的能量函数,以最小化能量目标函数为目标,控制参数曲线变形,具有最小能量的闭合曲线就是目标轮廓。

构造Snakes模型的目的是为了调和上层知识和底层图像特征这一对矛盾。无论是亮度、梯度、角点、纹理还是光流,所有的图像特征都是局部的。所谓局部性就是指图像上某一点的特征只取决于这一点所在的邻域,而与物体的形状无关。但是人们对物体的认识主要是来自于其外形轮廓。如何将两者有效地融合在一起正是Snakes模型的长处。Snakes模型的轮廓线承载了上层知识,而轮廓线与图像的匹配又融合了底层特征。这两项分别表示为Snakes模型中能量函数的内部力和图像力。

模型的形变受到同时作用在模型上的许多不同的力所控制,每一种力所产生一部分能量,这部分能量表示为活动轮廓模型的能量函数的一个独立的能量项。

Snake模型首先需要在感兴趣区域的附近给出一条初始曲线,接下来最小化能量泛函,让曲线在图像中发生变形并不断逼近目标轮廓。

Kass等提出的原始Snakes模型由一组控制点:v(s)=[x(s), y(s)] s∈[0, 1] 组成,这些点首尾以直线相连构成轮廓线。其中x(s)和y(s)分别表示每个控制点在图像中的坐标位置。 s 是以傅立叶变换形式描述边界的自变量。在Snakes的控制点上定义能量函数(反映能量与轮廓之间的关系):
其中第1项称为弹性能量是v的一阶导数的模,第2项称为弯曲能量,是v的二阶导数的模,第3项是外部能量(外部力),在基本Snakes模型中一般只取控制点或连线所在位置的图像局部特征例如梯度:
也称图像力。(当轮廓C靠近目标图像边缘,那么C的灰度的梯度将会增大,那么上式的能量最小,由曲线演变公式知道该点的速度将变为0,也就是停止运动了。这样,C就停在图像的边缘位置了,也就完成了分割。那么这个的前提就是目标在图像中的边缘比较明显了,否则很容易就越过边缘了。)

弹性能量和弯曲能量合称内部能量(内部力),用于控制轮廓线的弹性形变,起到保持轮廓连续性和平滑性的作用。而第三项代表外部能量,也被称为图像能量,表示变形曲线与图像局部特征吻合的情况。内部能量仅仅跟snake的形状有关,而跟图像数据无关。而外部能量仅仅跟图像数据有关。在某一点的α和β的值决定曲线可以在这一点伸展和弯曲的程度。

最终对图像的分割转化为求解能量函数Etotal(v)极小化(最小化轮廓的能量)。在能量函数极小化过程中,弹性能量迅速把轮廓线压缩成一个光滑的圆,弯曲能量驱使轮廓线成为光滑曲线或直线,而图像力则使轮廓线向图像的高梯度位置靠拢。基本Snakes模型就是在这3个力的联合作用下工作的。

因为图像上的点都是离散的,所以我们用来优化能量函数的算法都必须在离散域里定义。所以求解能量函数Etotal(v)极小化是一个典型的变分问题(微分运算中,自变量一般是坐标等变量,因变量是函数;变分运算中,自变量是函数,因变量是函数的函数,即数学上所谓的泛函。对泛函求极值的问题,数学上称之为变分法)。

在离散化条件(数字图像)下,由欧拉方程可知最终问题的答案等价于求解一组差分方程:(欧拉方程是泛函极值条件的微分表达式,求解泛函的欧拉方程,即可得到使泛函取极值的驻函数,将变分问题转化为微分问题。)


记外部力 F = −∇ P, Kass等将上式离散化后,对x(s)和y(s)分别构造两个五对角阵的线性方程组,通过迭代计算进行求解。在实际应用中一般先在物体周围手动点出控制点作为Snakes模型的起始位置,然后对能量函数迭代求解。


Reference:


李天庆等,Snake模型综述,计算机工程,2005,第31卷 第9期


代码实现

% =========================================================================  
%                   Snakes:Active Contour Models  
% =========================================================================  
% By gujinjin 2012/12/10-12  Sunny  
% 基于KASS等的论文思想  
% 参考文献:  
% [1] KASS etc.Snakes:Active Contour Models  
% [2] CSDN 博客 - Author:乐不思蜀Tone  
% [3] Ritwik Kumar(Harvard University),D.Kroon(Twente University)的程序  
% [4] 《数学建模与数学实验》  
% ------  
clc;clf;clear all;  
  
% =========================================================================  
%                      获取手动取点坐标  
% =========================================================================  
% 读取显示图像  
%I = imread('Coronary_MPR.jpg');  
I =  imread('test3.png');  
% 转化为双精度型  
%I = im2double(I);   
% 若为彩色,转化为灰度  
if(size(I,3)==3), I=rgb2gray(I); end  
figure(1),imshow(I);  
%---------------------------  
%        高斯滤波  
%---------------------------  
sigma=1.0;  
% 创建特定形式的二维高斯滤波器H  
H = fspecial('gaussian',ceil(3*sigma), sigma);  
% 对图像进行高斯滤波,返回和I等大小矩阵  
Igs = filter2(H,I,'same');  
%---------------------------  
%      获取Snake的点坐标  
%---------------------------  
figure(2),imshow(Igs);  
x=[];y=[];c=1;N=100; %定义取点个数c,上限N  
% 获取User手动取点的坐标  
% [x,y]=getpts  
while c<N  
    [xi,yi,button]=ginput(1);  
    % 获取坐标向量  
    x=[x xi];  
    y=[y yi];  
    hold on  
    % text(xi,yi,'o','FontSize',10,'Color','red');  
    plot(xi,yi,'ro');  
    % 若为右击,则停止循环  
    if(button==3), break; end  
    c=c+1;  
end  
  
% 将第一个点复制到矩阵最后,构成Snake环  
xy = [x;y];  
c=c+1;  
xy(:,c)=xy(:,1);  
% 样条曲线差值  
t=1:c;  
ts = 1:0.1:c;  
xys = spline(t,xy,ts);  
xs = xys(1,:);  
ys = xys(2,:);  
% 样条差值效果  
hold on  
temp=plot(x(1),y(1),'ro',xs,ys,'b.');  
legend(temp,'原点','插值点');  
  
% =========================================================================  
%                     Snakes算法实现部分  
% =========================================================================  
NIter =1000; % 迭代次数  
alpha=0.2; beta=0.2; gamma = 1; kappa = 0.1;  
wl = 0; we=0.4; wt=0;  
[row col] = size(Igs);  
  
% 图像力-线函数  
Eline = Igs;  
% 图像力-边函数  
[gx,gy]=gradient(Igs);  
Eedge = -1*sqrt((gx.*gx+gy.*gy));  
% 图像力-终点函数  
% 卷积是为了求解偏导数,而离散点的偏导即差分求解  
m1 = [-1 1];   
m2 = [-1;1];  
m3 = [1 -2 1];   
m4 = [1;-2;1];  
m5 = [1 -1;-1 1];  
cx = conv2(Igs,m1,'same');  
cy = conv2(Igs,m2,'same');  
cxx = conv2(Igs,m3,'same');  
cyy = conv2(Igs,m4,'same');  
cxy = conv2(Igs,m5,'same');  
  
for i = 1:row  
    for j= 1:col  
        Eterm(i,j) = (cyy(i,j)*cx(i,j)*cx(i,j) -2 *cxy(i,j)*cx(i,j)*cy(i,j) + cxx(i,j)*cy(i,j)*cy(i,j))/((1+cx(i,j)*cx(i,j) + cy(i,j)*cy(i,j))^1.5);  
    end  
end  
  
%figure(3),imshow(Eterm);  
%figure(4),imshow(abs(Eedge));  
% 外部力 Eext = Eimage + Econ  
Eext = wl*Eline + we*Eedge + wt*Eterm;  
% 计算梯度  
[fx,fy]=gradient(Eext);  
  
xs=xs';  
ys=ys';  
[m n] = size(xs);  
[mm nn] = size(fx);  
  
% 计算五对角状矩阵  
% 附录: 公式(14) b(i)表示vi系数(i=i-2 到 i+2)  
b(1)=beta;  
b(2)=-(alpha + 4*beta);  
b(3)=(2*alpha + 6 *beta);  
b(4)=b(2);  
b(5)=b(1);  
  
A=b(1)*circshift(eye(m),2);  
A=A+b(2)*circshift(eye(m),1);  
A=A+b(3)*circshift(eye(m),0);  
A=A+b(4)*circshift(eye(m),-1);  
A=A+b(5)*circshift(eye(m),-2);  
  
% 计算矩阵的逆  
[L U] = lu(A + gamma.* eye(m));  
Ainv = inv(U) * inv(L);   
  
% =========================================================================  
%                      画图部分  
% =========================================================================  
%text(10,10,'+','FontName','Time','FontSize',20,'Color','red');  
% 迭代计算  
figure(3)  
for i=1:NIter;  
    ssx = gamma*xs - kappa*interp2(fx,xs,ys);  
    ssy = gamma*ys - kappa*interp2(fy,xs,ys);  
   
    % 计算snake的新位置  
    xs = Ainv * ssx;  
    ys = Ainv * ssy;  
      
    % 显示snake的新位置  
    imshow(I);   
    hold on;  
    plot([xs; xs(1)], [ys; ys(1)], 'r-');  
    hold off;  
    pause(0.001)      
end  

运行效果



snake 模型

阅读数 4437

没有更多推荐了,返回首页