水平集LevelSet的使用問題
阿新 • • 發佈：2019-02-20
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); }