基本概念
OpenCV之前有个cvSnakeImage函数实现了Snakes活动轮廓模型，但是cvSnakeImage函数在opencv2中已经被去掉。由于没有无法得知其依据的论文，这里根据源代码简单介绍它的原理。遍历每个轮廓的点，计算该点邻域内每个点的能量E，将能量最小的点代替当前点。E= alpha * Econt + beta * Ecurv + gamma * Eimg。Econ表示轮廓点的距离信息，Ecurv表示轮廓点的曲率信息，Eimg表示像素或者像素梯度。
用OpenCV3.x重新实现后，发现貌似只对二值化图像分割有作用。
示例演示
用OpenCV3.x重新实现cvSnakeImage函数。完整工程代码。
void SnakeImage::segment(std::vector<cv::Point> &points)
{
drawCurve(points);
if (points.size() < 3)
return;
int i = 0, j = 0, k = 0;
//win_ is a domain, in which search minimization of energy for each point
//neighbors is the number of points in win_
int neighbors = win_.height * win_.width;
//the center of win_
int centerx = win_.width >> 1;
int centery = win_.height >> 1;
//energy terms of neighbors of every points
float *Econt = new float[neighbors];
float *Ecurv = new float[neighbors];
float *Eimg = new float[neighbors];
float *E = new float[neighbors];
//for two order difference
int tx = 0, ty = 0;
//for calculating distance
float averagedistance = 0, diffx = 0, diffy = 0;
//for calculating moving
int offsetx = 0, offsety = 0;
//temporary variables
float tmp = 0.0, energy = 0.0;
//if SNAKE_GRAD
cv::Mat dx, dy, src;
src_.convertTo(src, CV_16SC1);
Sobel(src, dx, src.depth(), 1, 0, 1);
Sobel(src, dy, src.depth(), 0, 1, 1);
// the number of moved points in each iteration
int movecount = 0;
bool converged = false;
int iteration = 0;
while (!converged)
{
movecount = 0; // the number of moved points
iteration++;
//the average of distance between points
averagedistance = 0;
for (i = 0; i < points.size(); i++)
{
if (i != 0)
{
diffx = points[i - 1].x - points[i].x;
diffy = points[i - 1].y - points[i].y;
}
else //n-1->0
{
diffx = points[points.size() - 1].x - points[i].x;
diffy = points[points.size() - 1].y - points[i].y;
}
averagedistance += std::sqrt(std::pow(diffx, 2) + std::pow(diffy, 2));
}
averagedistance /= points.size();
//search each point
for (i = 0; i < points.size(); i++)
{
//calculate Econt
float maxEcont = 0;
float minEcont = SNAKE_BIG;
//calculate reasonable search boundaries
int left = std::min(points[i].x, centerx);
int right = std::min(src_.cols - 1 - points[i].x, centerx);
int upper = std::min(points[i].y, centery);
int bottom = std::min(src_.rows - 1 - points[i].y, centery);
for (j = -upper; j <= bottom; j++)
{
for (k = -left; k <= right; k++)
{
if (i == 0)
{
diffx = points[points.size() - 1].x - (points[i].x + k);
diffy = points[points.size() - 1].y - (points[i].y + j);
}
else
{
diffx = points[i - 1].x - (points[i].x + k);
diffy = points[i - 1].y - (points[i].y + j);
}
energy = std::abs(averagedistance - std::sqrt(std::pow(diffx, 2) + std::pow(diffy, 2)));
Econt[(j + centery) * win_.width + k + centerx] = energy;
maxEcont = std::max(maxEcont, energy);
minEcont = std::min(minEcont, energy);
}
}
//if maxEcont == minEcont, set Econt all to zero
//else normalize Econt
tmp = maxEcont - minEcont;
tmp = (tmp == 0.0f) ? 0.0f : (1.0 / tmp);
for (k = 0; k < neighbors; k++)
{
Econt[k] = (Econt[k] - minEcont) * tmp;
}
//calculate Ecurv
float maxEcurv = 0.0f;
float minEcurv = SNAKE_BIG;
for (j = -upper; j <= bottom; j++)
{
for (k = -left; k <= right; k++)
{
if (i == 0)
{
tx = points[points.size() - 1].x - 2 * (points[i].x + k) + points[i + 1].x;
ty = points[points.size() - 1].y - 2 * (points[i].y + j) + points[i + 1].y;
}
else if (i == points.size() - 1)
{
tx = points[i - 1].x - 2 * (points[i].x + k) + points[0].x;
ty = points[i - 1].y - 2 * (points[i].y + j) + points[0].y;
}
else
{
tx = points[i - 1].x - 2 * (points[i].x + k) + points[i + 1].x;
ty = points[i - 1].y - 2 * (points[i].y + j) + points[i + 1].y;
}
energy = std::pow(tx, 2) + std::pow(ty, 2);
Ecurv[(j + centery) * win_.width + k + centerx] = energy;
maxEcurv = std::max(maxEcurv, energy);
minEcurv = std::min(minEcurv, energy);
}
}
// normalize as above
tmp = maxEcurv - minEcurv;
tmp = (tmp == 0.0f) ? 0.0f : (1.0 / tmp);
for (k = 0; k < neighbors; k++)
{
Ecurv[k] = (Ecurv[k] - minEcurv) * tmp;
}
//calculate Eimg
float maxEimg = 0.0f;
float minEimg = SNAKE_BIG;
for (j = -upper; j <= bottom; j++)
{
for (k = -left; k <= right; k++)
{
if (scheme_ == SNAKE_GRAD)
{
energy = std::pow(dx.at<short>(points[i].y + j, points[i].x + k), 2);
energy += std::pow(dy.at<short>(points[i].y + j, points[i].x + k), 2);
Eimg[(j + centery) * win_.width + k + centerx] = energy;
}
else
{
energy = src_.at<uchar>(points[i].y + j, points[i].x + k) ;
Eimg[(j + centery) * win_.width + k + centerx] = energy;
}
maxEimg = std::max(maxEimg, energy);
minEimg = std::min(minEimg, energy);
}
}
// normalize as above
tmp = (maxEimg - minEimg);
tmp = (tmp == 0.0f) ? 0.0f : (1.0 / tmp);
for (k = 0; k < neighbors; k++)
{
Eimg[k] = (minEimg - Eimg[k]) * tmp;
}
//calculate energy of neighbors
for (k = 0; k < neighbors; k++)
{
E[k] = alpha_ * Econt[k] + beta_ * Ecurv[k] + gamma_ * Eimg[k];
}
//find minimum energy of neighbors
float Emin = SNAKE_BIG;
for (j = -upper; j <= bottom; j++)
{
for (k = -left; k <= right; k++)
{
if (E[(j + centery) * win_.width + k + centerx] < Emin)
{
Emin = E[(j + centery) * win_.width + k + centerx];
offsetx = k;
offsety = j;
}
}
}
//if moved
if (offsetx || offsety)
{
points[i].x += offsetx;
points[i].y += offsety;
movecount++;
}
} // for each point
//judge convergence
converged = (movecount == 0);
if ((termcriteria_.type & CV_TERMCRIT_ITER) && (iteration >= termcriteria_.maxCount))
converged = 1;
if ((termcriteria_.type & CV_TERMCRIT_EPS) && (movecount <= termcriteria_.epsilon))
converged = 1;
drawCurve(points);
} //for each iteration
//free resource
delete[] Econt;
delete[] Ecurv;
delete[] Eimg;
delete[] E;
}
运行结果
相机
圆
来源：https://blog.csdn.net/webzhuce/article/details/100834086