几何活动轮廓模型——水平集分割：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 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();}
}

运行结果

参考资料

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

图像分割之水平集（Level Set）分割相关推荐

[图像处理]水平集(Level set)算法实现思路(简化)
[图像处理]水平集(Level set)算法实现思路(简化) 创建时间:2020年6月22日修改时间:2021年6月12日文章目录 [图像处理]水平集(Level set)算法实现思路(简化) 一 ...
水平集(level set)算法原理介绍
本篇文章,解释的是水平集算法最基础的原理. 1 水平集方法的解释有一个表面S,它与一个平面P相交,得到一个曲线C,这个C就是我们通过水平集得到的轮廓. 在图像分割中,表面S是随着由图像派生得到的 ...
基于水平集方法和G0模型的SAR图像分割
基于水平集方法和G0模型的SAR图像分割 Abstract(摘要) 这篇文章提出了一种分割SAR图像的方法,探索利用SAR数据中的统计特性将图像分区域.我们假设为SAR图像分割分配参数,并与水平集模型 ...
水平集（Level Set）的基本方法
水平集(Level Set)的基本方法水平集(Level Set)的基本方法-曲线演化的直观解释映射C(p), p\in [a,b] : R→R^2定义了一个平面的曲线,p是参数,对属于区间[a, ...
基于水平集的图像分割方法
一.引言借鉴一些流体中的重要思想, 1988年,Osher和Sethian首次提出了水平集算法[1],这是一种有效解决曲线演化问题的数值方法,并且计算稳定,适宜任意维数空间.随后,Osher等人对水 ...
matlab气管分割,一种基于区域生长法与水平集相融合的肺部CT图像的分割
摘要: 为将肺实质区域从含有背景.噪声的胸腔区域里分割出来,首先,应用传统的区域生长法初步定位肺部边界轮廓:其次,去除肺部边界噪声,采用自适应曲率阈值法修复肺部边界:最后,应用水平集法中的DRLSE模 ...
卷积神经网络结合水平集方法
水平集方法用于深度卷积网络水平集简介 CNN结合Level Set 疑惑水平集简介水平集(Level Set)方法是用于图像分割非常受欢迎的方法,其通过比目标维度高一维的水平集函数(LSF)的零 ...
图形学笔记（九）几何 ——几何表示方法（CSG、距离函数、水平集、点云、网格（obj格式））、贝塞尔曲线（面）
图形学笔记(八)着色2 -- 纹理映射.重心坐标.双线性插值.Mipmap.三线性插值.各向异性过滤.纹理的应用(环境贴图.法线贴图等) 图形学笔记(十)几何2 -- 曲面细分(Loop细分.Catm ...
荔枝hsv空间图像分割程序matlab,基于稀疏场水平集的荔枝图像分割算法_毛亮
第4期毛亮等:基于稀疏场水平集的荔枝图像分割算法349 从表中可知,本文算法在计算速度上接近标准模糊C 均值聚类算法,优于水平集方法.因此,与水平集方法相比,本文算法有着更好的分割性能和实时性. 4结 ...

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