水平集LevelSet的使用问题

1、水平集的是一个很好的分割算法，其使用比较多的是医学图像领域。其有一些有点也有一些缺点

其会计算图片的前景灰度、背景灰度。前景是灰度值大的目标，如果前景目标太小，则会导致其分割的不准确。其算法原理跟大津算法类似。其分割可以补偿一点缺失的边缘，这个特性很适合医学图像领域，因为医学图像很多噪点，边缘不明确的现象，还有一个优点是分割速度快。下面是一些例子：

分割成功并补偿了缺失边缘的：

分割失败：

灯光条太小，其分割成的是整张图

下面是水平集的原代码：

levelset.hpp

#include <vector>  
#include <opencv2/opencv.hpp> 
using namespace cv;
using namespace std;



class LevelSet
{
public:
	LevelSet();
	~LevelSet();

	//基本参数  
	int m_iterNum;      //迭代次数  
	float m_lambda1;    //全局项系数  
	float m_nu;     //长度约束系数ν  
	float m_mu;     //惩罚项系数μ  
	float m_timestep;   //演化步长δt  
	float m_epsilon;    //规则化参数ε  

	//过程数据  
	Mat m_mImage;       //源图像  

	int m_iCol;     //图像宽度  
	int m_iRow;     //图像高度  
	int m_depth;        //水平集数据深度  
	float m_FGValue;    //前景值  
	float m_BKValue;    //背景值  

	//初始化水平集  
	void initializePhi(Mat img,  //输入图像  
		int iterNum, //迭代次数  
		Rect boxPhi);//前景区域  
	void EVolution();   //演化  

	Mat m_mPhi;     //水平集：φ  
protected:
	Mat m_mDirac;       //狄拉克处理后水平集：δ（φ）  
	Mat m_mHeaviside;   //海氏函数处理后水平集：Н（φ）  
	Mat m_mCurv;        //水平集曲率κ=div(▽φ/|▽φ|)  
	Mat m_mK;       //惩罚项卷积核  
	Mat m_mPenalize;    //惩罚项中的▽<sup>2</sup>φ  

	void Dirac();       //狄拉克函数  
	void Heaviside();   //海氏函数  
	void Curvature();   //曲率  
	void BinaryFit();   //计算前景与背景值  
};

levelset.cpp

#include<iostream>
#include<opencv.hpp>
#include"levelset.hpp"
using namespace std;
using namespace cv;


LevelSet::LevelSet()
{
	m_iterNum = 300;
	m_lambda1 = 1;
	m_nu = 0.001 * 255 * 255;
	m_mu = 1.0;
	m_timestep = 0.1;
	m_epsilon = 1.0;
}


LevelSet::~LevelSet()
{
}

void LevelSet::initializePhi(Mat img, int iterNum, Rect boxPhi)
{
	//boxPhi是前景区域  
	m_iterNum = iterNum;
	cvtColor(img, m_mImage, CV_BGR2GRAY);

	m_iCol = img.cols;
	m_iRow = img.rows;
	m_depth = CV_32FC1;

	//显式分配内存  
	m_mPhi = Mat::zeros(m_iRow, m_iCol, m_depth);
	m_mDirac = Mat::zeros(m_iRow, m_iCol, m_depth);
	m_mHeaviside = Mat::zeros(m_iRow, m_iCol, m_depth);

	//初始化惩罚性卷积核  
	m_mK = (Mat_<float>(3, 3) << 0.5, 1, 0.5,
		1, -6, 1,
		0.5, 1, 0.5);

	int c = 2;
	for (int i = 0; i < m_iRow; i++)
	{
		for (int j = 0; j < m_iCol; j++)
		{
			if (i<boxPhi.y || i>boxPhi.y + boxPhi.height || j<boxPhi.x || j>boxPhi.x + boxPhi.width)
			{
				m_mPhi.at<float>(i, j) = -c;
			}
			else
			{
				m_mPhi.at<float>(i, j) = c;
			}
		}
	}
}

void LevelSet::Dirac()
{
	//狄拉克函数  
	float k1 = m_epsilon / CV_PI;
	float k2 = m_epsilon*m_epsilon;
	for (int i = 0; i < m_iRow; i++)
	{
		float *prtDirac = &(m_mDirac.at<float>(i, 0));
		float *prtPhi = &(m_mPhi.at<float>(i, 0));

		for (int j = 0; j < m_iCol; j++)
		{
			float *prtPhi = &(m_mPhi.at<float>(i, 0));
			prtDirac[j] = k1 / (k2 + prtPhi[j] * prtPhi[j]);
		}
	}
}

void LevelSet::Heaviside()
{
	//海氏函数  
	float k3 = 2 / CV_PI;
	for (int i = 0; i < m_iRow; i++)
	{
		float *prtHeaviside = (float *)m_mHeaviside.ptr(i);
		float *prtPhi = (float *)m_mPhi.ptr(i);

		for (int j = 0; j < m_iCol; j++)
		{
			prtHeaviside[j] = 0.5 * (1 + k3 * atan(prtPhi[j] / m_epsilon));
		}
	}
}

void LevelSet::Curvature()
{
	//计算曲率  
	Mat dx, dy;
	Sobel(m_mPhi, dx, m_mPhi.depth(), 1, 0, 1);
	Sobel(m_mPhi, dy, m_mPhi.depth(), 0, 1, 1);

	for (int i = 0; i < m_iRow; i++)
	{
		float *prtdx = (float *)dx.ptr(i);
		float *prtdy = (float *)dy.ptr(i);
		for (int j = 0; j < m_iCol; j++)
		{
			float val = sqrtf(prtdx[j] * prtdx[j] + prtdy[j] * prtdy[j] + 1e-10);
			prtdx[j] = prtdx[j] / val;
			prtdy[j] = prtdy[j] / val;
		}
	}
	Mat ddx, ddy;
	Sobel(dx, ddy, m_mPhi.depth(), 0, 1, 1);
	Sobel(dy, ddx, m_mPhi.depth(), 1, 0, 1);
	m_mCurv = ddx + ddy;
}

void LevelSet::BinaryFit()
{
	//先计算海氏函数  
	Heaviside();

	//计算前景与背景灰度均值  
	float sumFG = 0;
	float sumBK = 0;
	float sumH = 0;
	//float sumFH = 0;  
	Mat temp = m_mHeaviside;
	Mat temp2 = m_mImage;
	float fHeaviside;
	float fFHeaviside;
	float fImgValue;
	for (int i = 1; i < m_iRow; i++)
	{
		float *prtHeaviside = &(m_mHeaviside.at<float>(i, 0));
		uchar *prtImgValue = &(m_mImage.at<uchar>(i, 0));
		for (int j = 1; j < m_iCol; j++)
		{
			fImgValue = prtImgValue[j];
			fHeaviside = prtHeaviside[j];
			fFHeaviside = 1 - fHeaviside;

			sumFG += fImgValue*fHeaviside;
			sumBK += fImgValue*fFHeaviside;
			sumH += fHeaviside;
		}
	}
	m_FGValue = sumFG / (sumH + 1e-10);         //前景灰度均值  
	m_BKValue = sumBK / (m_iRow*m_iCol - sumH + 1e-10); //背景灰度均值  
}
Mat showIMG;
void LevelSet::EVolution()
{
	float fCurv;
	float fDirac;
	float fPenalize;
	float fImgValue;

	for (int i = 0; i < m_iterNum; i++)
	{
		Dirac();
		Curvature();
		BinaryFit();
		filter2D(m_mPhi, m_mPenalize, m_depth, m_mK, Point(1, 1));//惩罚项的△φ  
		for (int i = 0; i < m_iRow; i++)
		{
			float *prtCurv = &(m_mCurv.at<float>(i, 0));
			float *prtDirac = &(m_mDirac.at<float>(i, 0));
			float *prtPenalize = &(m_mPenalize.at<float>(i, 0));
			uchar *prtImgValue = &(m_mImage.at<uchar>(i, 0));
			for (int j = 0; j < m_iCol; j++)
			{
				fCurv = prtCurv[j];
				fDirac = prtDirac[j];
				fPenalize = prtPenalize[j];
				fImgValue = prtImgValue[j];

				float lengthTerm = m_nu* fDirac * fCurv;                    //长度约束  
				float penalizeTerm = m_mu*(fPenalize - fCurv);                  //惩罚项  
				float areaTerm = fDirac * m_lambda1 *                       //全局项  
					(-((fImgValue - m_FGValue)*(fImgValue - m_FGValue))
					+ ((fImgValue - m_BKValue)*(fImgValue - m_BKValue)));

				m_mPhi.at<float>(i, j) = m_mPhi.at<float>(i, j) + m_timestep*(lengthTerm + penalizeTerm + areaTerm);
			}
		}

		//显示每一次演化的结果  


		//return showIMG;
	}
	cvtColor(m_mImage, showIMG, CV_GRAY2BGR);
	Mat Mask = m_mPhi >= 0;   //findContours的输入是二值图像  
	dilate(Mask, Mask, Mat(), Point(-1, -1), 3);
	erode(Mask, Mask, Mat(), Point(-1, -1), 3);
	vector<vector<Point> > contours;
	findContours(Mask,
		contours,// 轮廓点  
		RETR_EXTERNAL,// 只检测外轮廓  
		CHAIN_APPROX_NONE);// 提取轮廓所有点  
	drawContours(showIMG, contours, -1, Scalar(255, 0, 0), 2);
	namedWindow("Level Set后图像");
	imshow("Level Set后图像", showIMG);
	waitKey(1);
}
void main()
{
	Mat img = imread("0.jpg");
	imshow("原图", img);
	LevelSet ls;
	VideoCapture cap;
	cap.open(0);
	cap.set(CV_CAP_PROP_FRAME_WIDTH, 2048);
	cap.set(CV_CAP_PROP_FRAME_HEIGHT, 1536);
	for (;;)
	{
		cap >> img;
		resize(img, img, Size(img.cols / 8, img.rows / 8));
		imshow("原图", img);
		Rect rec(0, 0, img.cols, img.rows);
		ls.initializePhi(img, 30, rec);
		ls.EVolution();
	}
	//imshow("Level Set后图像", showIMG);
	waitKey(0);
}