【主动轮廓模型（⼀）】《Snakes ：ActiveContourModels 》算法原理与O 。。
。
⽂章⽬录
1 概述
主动轮廓模型（也称Active Contour Model、Snake）是Kass等⼈在1988年提出的，该算法将图像分割问题转换为求解能量泛函最⼩值的问题。主要思路是通过构造能量泛函，经过算法迭代，轮廓曲线由初始位置逐渐向使能量函数最⼩（或局部极⼩）的图像边缘逼近，最终分割出⽬标。
2 算法原理
⾸先需要⼈为地在图像上给出初始轮廓曲线，确切的说是⼀组⽤于控制曲线形状的控制点：，这些点收尾相连构成⼀个封闭的轮廓线。其中和分别表⽰每个控制点在图像中的坐标位置，是以傅⽴叶变换形式描述边界的⾃变量，也可以理解为弧长。则Snake曲线的能量函数表⽰为：
其中，为内部能量，为图像能量，为外部约束能量。其中图像能量和外部约束能量统称为外部能量，即：。
2.1 内部能量 内部能量由保证曲线连续性的⼀阶导和保证曲线平滑的⼆阶导组成，表⽰为：粽子蒸多久
通过调整权值和可以控制曲线的形状。例如将置为0可以让曲线最终出现拐⾓，即曲线⼆阶不连续。
2.2 图像能量 在低层次的计算机视觉中，需要能够将轮廓吸引到特定图像特征的能量函数。原始的Snake轮廓模型中提出了3种不同的能量函数，分别将snake轮廓吸引到线、边和末端。完整的图像能量 可以表⽰为这3个能量函数的权值组合，通过调整这3个权值，可以形成不同的轮廓形状。
（1）线函数最简单直接且有⽤的图像能量函数是图像本⾝，即图像本⾝的灰度。令：
v (s )=[x (s ),y (s )]s ∈[0,1]x (s )y (s )s E snake ∗
=E (v (s ))ds
∫01
snake =E (v (s ))+E (v (s ))+E (v (s ))ds
∫01
int image con =E (v (s ))+E (v (s ))ds
∫01
int ext (1)
E int E image E con E =
ext E +image E con E int
E int E =int (α(s )∣v (s )∣+21
s 2β(s )∣v (s )∣)
ss 2(2)
α(s )β(s )β(s )E image
E image E =image ωE +line line ωE +edge edge ωE term term (3)
E line
控制的正负号可以控制轮廓被吸引到较暗的线或是较亮的线，也就是使轮廓试图靠近轮廓的最暗或最亮处。
然⽽，如果snake轮廓的⼀部分到达了⼀个低能量的图像特征位置，这个样条项将推动snake轮廓临近的部分朝着这个特征可能的延续⽅向移动，这会使其在⼀个最优的局部最⼩位置引⼊⼀个较⼤的能量。⼀种解决⽅案是允许snake轮廓和模糊能量函数平衡，然后慢慢降低模糊程度。
Marr和Hildreth在他们的论⽂⾥证明了“灰度的突然变化会在⼀阶导数中引起波峰或波⾕，或在⼆阶导数中等效地引起零交点”。为了显⽰图像尺度空间连续性和Marr-Hildreth边缘检测理论的关系，snake模型中采⽤了模糊的边能量函数：
其中是标准差为的⾼斯函数，该函数的最⼩值位于在Marr-Hildreth理论中被定义的的零交点处。在能量函数中加⼊这⼀项意味着snake轮廓在被吸引到零交点的同时，仍然受到它⾃⼰的平滑限制。
（2）边函数在图像上找边缘可以通过梯度来实现，令：
在某个点上，梯度越⼤上式的能量越⼩，则snake轮廓将被吸引到梯度较⼤的区域。
（3）末端函数为了找到轮廓的终⽌位置，作者将平滑过图像中等⾼线的曲率加⼊到能量函数中。令是模糊后的图像，是梯度⾓，、是沿着和垂直梯度⽅向的单位向量。则中等⾼线的曲率可以表⽰
为：
通过组合和，我们能够创建⼀个被边和末端吸引的snake曲线。
2.3 外部能量 外部能量来⾃于外部的约束⼒。在论⽂原⽂中，作者给了⼀个外部约束的例⼦， 包括外部的固定点、连接两条snake曲线的锚点或⿏标拖动的点。例如，为了在和点之间创建⼀个连接的弹性⼒，就可以将添加到外部能量中。
E =line I (x ,y )
(4)
ωline E =line −(G ∗σ∇I )22
(6)
G σσG ∗σ∇I 2E edge
E =edge −∣∇I (x ,y )∣2
(5)
E term
C (x ,y )=G (x ,y )∗σI (x ,y )θ=
密切关注tan ()−1C x C y
n =(cosθ,sinθ)n =⊥(−sinθ,cosθ)C (x ,y )E term =∂n ⊥
∂θ
=∂C /∂n ∂C /∂n 2⊥
2
=
(C +C )x 2y 2
3/2
C C −2C C C +C C yy x 2
xy x y xx y
2(6)
E edge E term E con
E con x 1x 2−k (x −1x )22E con
一棵大白菜
3 模型求解
由欧拉⽅程，求解能量的最⼩值（局部极⼩值），公式（1）的导数必须满⾜：
由，将上式改写成x和y两个⽅向有：
利⽤顺序连接的锚点的坐标差分来近似导数：
将（9）带⼊（8），同时令有：
令：
E snake ∗
αv (s )−ss βv (s )−ssss ∇E (v (s ))−image ∇E (v (s ))=con 0
(7)
v (s )=[x (s ),y (s )]{αx (s )−βx (s )−=0
ss ssss ∂x ∂E ext
αy (s )−βy (s )−=0ss
ssss
∂y ∂E ext
(8)
⎩
⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧x (s )=x (s +1)+x (s −1)−2x (s )ss x (s )=ssss (x (s +2)+x (s )−2x (s +1))+(x (s )+x (s −2)−2x (s −1))−2(x (s +1)+x (s −1)−2x (s ))(9)
f =
x ,f =
∂x ∂E ext
y
∂y ∂E ext
{βx (s −2)−(α+4β)x (s −1)+(2α+6β)x (s )−(α+4β)x (s +1)+βx (s +2)+f =0
x βy (s −2)−(α+4β)y (s −1)+(2α+6β)y (s )−(α+4β)y (s +1)+βy (s +2)+f =0
y
(10)
则将（10）写成矩阵形式为：
其中A为五对⾓带状矩阵：
利⽤梯度下降法，在第t次迭代有：
其中是迭代步长，公式（14）的解为：
由于为五对⾓带状矩阵，因此也是⼀个五对⾓带状矩阵，原⽂中作者⽤LU分解来求其逆矩阵。
4 算法实现（OpenCV3）
⽹上关于snake算法的实现多是matlab代码。OpenCV2中可⽤cvSnakeImage来实现，但这个函数在OpenCV3中已经被删除。本⽂将在OpenCV3.14中实现snake算法代码。
cv ::Mat Interate ( cv ::Mat image , cv ::Mat xs , cv ::Mat ys , double alpha , double beta , double gamma , double kappa , double wl ,
⎩
⎪⎨⎪⎧a =2α+6βb =−(α+4β)c =β
(11)
{Ax +f (x ,y )=0
x Ay +f (x ,y )=0
y
(12)
A =,x =⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡a b c ⋮⋯c b八万里
b a b ⋮
c ⋯c c b a ⋮b c ⋯⋯c b ⋮a b c c ⋯c ⋮b a b b c ⋯⋮c b a ⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤,f =⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡x 1x 2x 3⋮x n −1x n x 1⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤x ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡f x 1f x 2f x 3⋮f x n −1f x n f x 1⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥
⎤(13)
{Ax +f (x ,y )=−γ(x −x )
t x t −1t −1t t −1Ay +f (y ,y )=−γ(y −y )
t y
t −1t −1t t −1(14)
γ{x =(A +γI )(x −f (x ,y ))
t −1t −1x t −1t −1y =(A +γI )(y −f (x ,y ))
t −1t −1y
t −1t −1(15)
A A +γI
double wl,
double we,
double wt,
int iterations
)
{
// 相关参数
int N = iterations;
cv::Mat smth = image.clone();
// 图像⼤⼩
qDebug()<<"Calculating size of image";
cv::Size size = image.size();
int row = size.height;
int col = size.width;
// 计算外部⼒（图像⼒）
qDebug()<<"Computing external forces";
cv::Mat E_line = smth.clone();// E_line is simply the image intensities
cv::Mat gradx, grady;
cv::Sobel(smth, gradx, smth.depth(),1,0,1,1,0, cv::BORDER_CONSTANT);
cv::Sobel(smth, grady, smth.depth(),0,1,1,1,0, cv::BORDER_CONSTANT);
qDebug()<<"Computing gradx and grady";
cv::Mat E_edge(row, col, CV_32FC1);
for(int i =0; i < ws; i++)
人皆有之{
for(int j =0; j < ls; j++)
{
float v_gradx = gradx.at<float>(i, j);
float v_grady = grady.at<float>(i, j);
E_edge.at<float>(i, j)=-1* std::sqrt(v_gradx * v_gradx + v_grady * v_grady);// E_edge is measured by gradient in the image
}
}
// 导数mask
qDebug()<<"masks for taking various derivatives";
cv::Mat m1 =(cv::Mat_<float>(1,2)<<-1,1);
cv::Mat m2 =(cv::Mat_<float>(2,1)<<-1,1);
cv::Mat m3 =(cv::Mat_<float>(1,3)<<1,-2,1);
cv::Mat m4 =(cv::Mat_<float>(3,1)<<1,-2,1);
cv::Mat m5 =(cv::Mat_<float>(2,2)<<1,-1,-1,1);
cv::Mat cx, cy, cxx, cyy, cxy;
filter2D(smth, cx,-1, m1);
filter2D(smth, cy,-1, m2);
宁的成语filter2D(smth, cxx,-1, m3);
filter2D(smth, cyy,-1, m4);
filter2D(smth, cxy,-1, m5);
// 计算 E_term
cv::Mat E_term(row, col, CV_32FC1);
for(int i =0; i < row; i++)
{
for(int j =0; j < col; j++)
祝学生学业有成的句子{
摊饼int v_cx = cx.at<float>(i, j);
int v_cy = cy.at<float>(i, j);
int v_cxx = cxx.at<float>(i, j);
int v_cyy = cyy.at<float>(i, j);
int v_cxy = cxy.at<float>(i, j);
E_term.at<float>(i, j)=(v_cyy*v_cx*v_cx -2* v_cxy*v_cx*v_cy + v_cxx * v_cy*v_cy)/(std::pow((1+ v_cx * v_cx + v_cy * v_cy),1.5)); }
