水平集方法

水平集是跟踪轮廓和表面运动的一种数字化方法，它不直接对轮廓进行操作，而是将轮廓设置成一个高维函数的零水平集。这个高维函数叫做水平集函数。然后对该水平集函数进行微分，通过从输出中提取零水平集来得到运动的轮廓。使用水平集的主要优点是可以对任意复杂的结构进行模式化和拓扑变换。
水平集（Level Set）方法是由Sethian和Osher于1988年提出。简单来说，水平集方法把低维的一些计算上升到更高一维，把N维的描述看成是N+1维的一个水平。比如，一个二维平面的圆，如x^2+y2=1，可以看成三维中二元函数f(x,y) = x^2+y2的1的水平集。如下图
在这里插入图片描述
专业描述：水平集方法将平面闭合曲线隐含地表达为连续函数曲面

几何活动轮廓模型

这里具体介绍2001年的文献《Active Contours Without Edges》，即Chan-Vese模型。作者提出了基于Mumford–Shah分割技术和水平集方法的活动轮廓模型，该模型不依赖边缘函数，即图像梯度。能量函数有fitting项和一些正则项组成。
在这里插入图片描述
作者把上面关于曲线C的能量函数写成关于水平集函数ϕ的能量函数，然后求解。

其中Heaviside函数和Delata函数，如下

最后得到水平集函数的迭代公式：

示例演示

下面基于OpenCV实现该算法。完整工程代码。
头文件

/**********************************************************************Copyright (c) Mr.Bin. All rights reserved.
For more information visit: http://blog.csdn.net/webzhuce**********************************************************************/
/**
* @brief   This is a filter that implements  Tony F. Chan. Active 
*          Contours Without Edges[J].2001
*/
#pragma once
#include "opencv2/opencv.hpp"class LevelSet
{
public:LevelSet(const cv::Mat &src);//initialize level setvoid initializePhi(cv::Point2f center, float radius);void evolving();	private:int iterationnum_ = 100;float lambda1_ = 1.0f;	 //weight for c1 fitting termfloat lambda2_ = 1.0f;   //weight for c2 fitting termfloat mu_ = 0.1 * 255 * 255;	//μ: weight for length termfloat nu_ = 0.0;  //ν: weight for area term, default value 0 float timestep_ = 0.1; //time step: δt//ε: parameter for computing smooth Heaviside and dirac functionfloat epsilon_ = 1.0;	float c1_;	//average(u0) inside level setfloat c2_;	//average(u0) outside level setcv::Mat phi_;		//level set: φcv::Mat dirac_;		//dirac level set: δ(φ)cv::Mat heaviside_;	//heaviside：Н(φ)cv::Mat curv_;		cv::Mat src_;cv::Mat image_; //for showingvoid gradient(const cv::Mat &src, cv::Mat &gradx, cv::Mat &grady);//Dirac functionvoid dirac();//Heaviside functionvoid heaviside();void calculateCurvature();//calculate c1 and c2void calculatC();//show evolvingvoid showEvolving();
};

源文件

/**********************************************************************Copyright (c) Mr.Bin. All rights reserved.
For more information visit: http://blog.csdn.net/webzhuce**********************************************************************/
#include "levelset.h"using namespace cv;
using namespace std;LevelSet::LevelSet(const cv::Mat &src)
{if (src.channels() == 3){cv::cvtColor(src, src_, CV_BGR2GRAY);src.copyTo(image_);}else{src.copyTo(src_);cv::cvtColor(src, image_, CV_GRAY2BGR);}src_.convertTo(src_, CV_32FC1);phi_ = Mat::zeros(src_.size(), CV_32FC1);dirac_ = Mat::zeros(src_.size(), CV_32FC1);heaviside_ = Mat::zeros(src_.size(), CV_32FC1);
}void LevelSet::initializePhi(cv::Point2f center, float radius)
{const float c = 2.0f;float value = 0.0;for (int i = 0; i < src_.rows; i++){for (int j = 0; j < src_.cols; j++){value = -sqrt(pow((j - center.x), 2) + pow((i - center.y), 2)) + radius;if (abs(value) < 0.1) //phi_.at<float>(i, j) = 0;else if (value >=0.1) //insidephi_.at<float>(i, j) = c;elsephi_.at<float>(i, j) = -c;}}
}void LevelSet::gradient(const cv::Mat &src, cv::Mat &gradx, cv::Mat &grady)
{if (src.rows < 2 || src.cols < 2)return;src.copyTo(gradx);src.copyTo(grady);for (int i = 0; i < src.rows; i++){for (int j = 0; j < src.cols; j++){if (j == 0)gradx.at<float>(i, j) = src.at<float>(i, j + 1) - src.at<float>(i, j);else if (j == src.cols - 1)gradx.at<float>(i, j) = src.at<float>(i, j) - src.at<float>(i, j - 1);elsegradx.at<float>(i, j) = (src.at<float>(i, j + 1) - src.at<float>(i, j - 1)) / 2.0;}}for (int j = 0; j < src.cols; j++){for (int i = 0; i < src.rows; i++){if (i == 0)grady.at<float>(i, j) = src.at<float>(i + 1, j) - src.at<float>(i, j);else if (i == src.rows - 1)grady.at<float>(i, j) = src.at<float>(i, j) - src.at<float>(i - 1, j);elsegrady.at<float>(i, j) = (src.at<float>(i + 1, j) - src.at<float>(i - 1, j)) / 2.0;}}}void LevelSet::dirac()
{const float k1 = epsilon_ / CV_PI;const float k2 = epsilon_*epsilon_;for (int i = 0; i < src_.rows; i++){for (int j = 0; j < src_.cols; j++){dirac_.at<float>(i, j) = k1 / (k2 + pow(phi_.at<float>(i, j),2));}}
}void LevelSet::heaviside()
{const float k3 = 2 / CV_PI;for (int i = 0; i < src_.rows; i++){for (int j = 0; j < src_.cols; j++){heaviside_.at<float>(i, j) = 0.5 * (1 + k3 * atan(phi_.at<float>(i, j) / epsilon_));}}
}void LevelSet::calculateCurvature()
{Mat dx, dy;gradient(src_, dx, dy);Mat norm;pow(dx.mul(dx) + dy.mul(dy), 0.5, norm);Mat dxx, dxy, dyx, dyy;gradient(dx / norm, dxx, dxy);gradient(dy / norm, dyx, dyy);curv_ = dxx + dyy;
}void LevelSet::calculatC()
{c1_ = 0.0f;c2_ = 0.0f;float sum1 = 0.0f;float h1 = 0.0f;float sum2 = 0.0f;float h2 = 0.0f;int outsidenum = 0;float value = 0.0f;float h = 0.0f;for (int i = 0; i < src_.rows; i++){for (int j = 0; j < src_.cols; j++){value = src_.at<float>(i, j);h = heaviside_.at<float>(i, j);h1 += h;sum1 = h *  value;h2 += (1 - h);sum2 += (1 - h) * value;}}c1_ = sum1 / (h1 + 1e-10);c2_ = sum2 / (h2 + 1e-10);
}void LevelSet::showEvolving()
{Mat image = image_.clone();Mat mask = phi_ >= 0;cv::dilate(mask, mask, Mat(), Point(-1, -1), 3);cv::erode(mask, mask, Mat(), Point(-1, -1), 3);vector<vector<Point> > contours;findContours(mask, contours, RETR_EXTERNAL, CHAIN_APPROX_NONE);drawContours(image, contours, -1, CV_RGB(0, 255,0), 2);imshow("Evolving", image);waitKey(0);
}void LevelSet::evolving()
{showEvolving();//iterationfor (int i = 0; i < iterationnum_; i++){heaviside();dirac();calculateCurvature();calculatC();//update phifor (int i = 0; i < src_.rows; i++){for (int j = 0; j < src_.cols; j++){float curv = curv_.at<float>(i, j);float dirac = dirac_.at<float>(i, j);float u0 = src_.at<float>(i, j);;float lengthTerm = mu_* dirac * curv;float areamterm = nu_ * dirac;float fittingterm = dirac * (-lambda1_ * pow(u0 - c1_, 2) + lambda2_ * pow(u0 - c2_, 2));float term = lengthTerm + areamterm + fittingterm;phi_.at<float>(i, j) = phi_.at<float>(i, j) + timestep_ * term;}}//just for showingshowEvolving();}
}