%?基本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為循環體的