基本Snake活动轮廓模型

% 基本Snake活动轮廓模型

I=imread（‘H:/MRI1/12.jpg‘）;   % 读入的图片应为uint8类型二维的灰度图
snake（I）;			    % 对图像I求其中需要分割物体的snake边界
   
function snake（I）
% Snake主体部分

alpha=0.5; beta=0;      % 连续参数alpha=0.5；平滑参数beta=0；步长为1
[xy]=DrawLine（I）;	% 在图像I上手动画线，得到初始轮廓线

a=2*alpha+6*beta; b=-（alpha+4*beta）; c=beta; 
J=[c b a b c]; h=max（size（x））;
A=diagCyclMat（hJ）;     % 求取设定参数下的五对角循环矩阵

II=eye（h）; [m~]=size（I）;        % 初始化
I=double（I）;

I1=-ff（I）;                       % 高斯势能I1
[I2xI2y]=NGradient（I1）;         % I1的负梯度I2
T=max（max（abs（I2x（:）））max（abs（I2y（:））））;
I2x=I2x/T; I2y=I2y/T;            % 梯度归一化
fx=-1*I2x; fy=-1*I2y;            % f为图像I的高斯势能的梯度


for t=1:10000                      % 迭代，未计算迭代终点   
    ffx=fx（m*（uint16（x）-1）+uint16（y））;
    ffy=fy（m*（uint16（x）-1）+uint16（y））;
    x=（（II/（A+II））*（x‘-ffx‘））‘;
    y=（（II/（A+II））*（y‘-ffy‘））‘;
end

I=uint8（I）; imshow（I）;  hold on
plot（xy‘Color‘‘White‘）         % 显示最终Snake轮廓线
end

function I1=ff（I）
%求取I的边缘函数（负高斯势能）

%5阶Standard Deviation=3的高斯滤波，sobel梯度
h=fspecial（‘gaussian‘53）; w1=fspecial（‘sobel‘）; w2=w1‘;

Is=imfilter（double（I）h‘conv‘‘replicate‘）;
I1=imfilter（Isw1‘replicate‘）.^2+imfilter（Isw2‘replicate‘）.^2;

end

function [I2xI2y]=NGradient（I）
%求取I的负梯度
%sobel梯度
w1=fspecial（‘sobel‘）; w2=w1‘;
I=double（I）;
I2y=imfilter（Iw1‘replicate‘）;
I2x=imfilter（Iw2‘replicate‘）;
end

function A=diagCyclMat（nJ）
% A = diagonal cycle（J） matrix.
%生成一个以向量J为循环体的