图像分割之水平集（Level Set）分割

几何活动轮廓模型——水平集分割：Active Contours Without Edges

水平集方法

        水平集是跟踪轮廓和表面运动的一种数字化方法，它不直接对轮廓进行操作，而是将轮廓设置成一个高维函数的零水平集。这个高维函数叫做水平集函数。然后对该水平集函数进行微分，通过从输出中提取零水平集来得到运动的轮廓。使用水平集的主要优点是可以对任意复杂的结构进行模式化和拓扑变换。
        水平集（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 set
	void initializePhi(cv::Point2f center, float radius);
	void evolving();	

private:
	int iterationnum_ = 100;
	float lambda1_ = 1.0f;	 //weight for c1 fitting term
	float lambda2_ = 1.0f;   //weight for c2 fitting term
	float mu_ = 0.1 * 255 * 255;	//μ: weight for length term
	float nu_ = 0.0;  //ν: weight for area term, default value 0 
	float timestep_ = 0.1; //time step: δt
	//ε: parameter for computing smooth Heaviside and dirac function
	float epsilon_ = 1.0;	
	float c1_;	//average(u0) inside level set
	float c2_;	//average(u0) outside level set
	cv::Mat phi_;		//level set: φ
	cv::Mat dirac_;		//dirac level set: δ(φ)
	cv::Mat heaviside_;	//heaviside：Н(φ)
	cv::Mat curv_;		
	cv::Mat src_;
	cv::Mat image_; //for showing
	
	void gradient(const cv::Mat &src, cv::Mat &gradx, cv::Mat &grady);
	//Dirac function
	void dirac();
	//Heaviside function
	void heaviside();
	void calculateCurvature();
	//calculate c1 and c2
	void calculatC();
	//show evolving
	void 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) //inside
				phi_.at<float>(i, j) = c;
			else
				phi_.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);
			else
				gradx.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);
			else
				grady.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();
	//iteration
	for (int i = 0; i < iterationnum_; i++)
	{
		heaviside();
		dirac();
		calculateCurvature();
		calculatC();

		//update phi
		for (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 showing
		showEvolving();
	}
}

运行结果

参考资料

• Level set method: Explanation
• 计算机视觉之图像分割——水平集方法_ACWE2001
• 水平集图像分割序列——CV模型