本文介绍了一种基于水平集的CV模型实现图像分割的方法，详细解析了从XML文件读取参数，初始化零水平集框，到使用狄拉克、海氏函数处理水平集，计算曲率，并进行演化的过程。
本项目使用水平集的CV模型
应用场景：分割图像
main.cpp
#include<iostream>
#include<opencv2/opencv.hpp>
#include"levelset.h"
//打开xml文件需要加载的头文件
#include "../tinyxml/tinyxml.h"
#include "../tinyxml/tinystr.h"
#include<string>
#include<vector>
using namespace std;
using namespace cv;
//定义初始化零水平集框
Rect boxPhi;
//保存解读xml获取rect
struct BoxSize
{
int xMin;
int yMin;
int xMax;
int yMax;
};
//xml文件解读
bool ReadParaXml(string m_strXmlPath, vector<BoxSize>& vecNode)
{
BoxSize* pNode = new BoxSize;
//读取xml文件中的参数值
TiXmlDocument* Document = new TiXmlDocument();
if (!Document->LoadFile(m_strXmlPath.c_str()))
{
cout << "无法加载xml文件！" << endl;
cin.get();
return false;
}
TiXmlElement* RootElement = Document->RootElement(); //根目录
TiXmlElement* NextElement = RootElement->FirstChildElement(); //根目录下的第一个节点层
while (NextElement != NULL) //判断有没有读完
{
if (NextElement->ValueTStr() == "object") //读到object节点
{
TiXmlElement* BoxElement = NextElement->FirstChildElement();
while (BoxElement->ValueTStr() != "bndbox") //读到box节点
{
BoxElement = BoxElement->NextSiblingElement();
}
//索引到xmin节点
TiXmlElement* xminElemeng = BoxElement->FirstChildElement();
{
//分别读取四个数值
pNode->xMin = atof(xminElemeng->GetText());
TiXmlElement* yminElemeng = xminElemeng->NextSiblingElement();
pNode->yMin = atof(yminElemeng->GetText());
TiXmlElement* xmaxElemeng = yminElemeng->NextSiblingElement();
pNode->xMax = atof(xmaxElemeng->GetText());
TiXmlElement* ymaxElemeng = xmaxElemeng->NextSiblingElement();
pNode->yMax = atof(ymaxElemeng->GetText());
//加入到向量中
vecNode.push_back(*pNode);
}
}
NextElement = NextElement->NextSiblingElement();
}
//释放内存
delete pNode;
delete Document;
cout << "完成xml的读取" << endl;
return true;
}
int main()
{
String xml_paths = "E:\\哈工智能\\fenge\\level_set\\tinyxml\\doc\\*.xml";
std::vector<cv::String> xml_files;
glob(xml_paths, xml_files);
String Img_paths = "C:\\Users\\tudejiang\\Desktop\\超声图像VOC格式标准数据集201911\\VOC2007\\JPEGImages\\*.jpg";
vector<String> image_files;
glob(Img_paths, image_files);
for (int i = 0; i < image_files.size(); i++)
{
vector<BoxSize> vecNode;
ReadParaXml(xml_files[i], vecNode);
boxPhi.x = vecNode[0].xMin+10 ;
boxPhi.y = vecNode[0].yMin+10 ;
boxPhi.width = vecNode[0].xMax - vecNode[0].xMin-20 ;
boxPhi.height = vecNode[0].yMax - vecNode[0].yMin-10;
cout << image_files[i]<< endl;
cout << xml_files[i] << endl;
Mat src = imread(image_files[i], 0);
//imshow("src", src);
LevelSet mylevelset;
mylevelset.initializePhi(src, 100, boxPhi);
double time0 = static_cast<double>(getTickCount());
mylevelset.EVolution(image_files[i]);
time0 = ((double)getTickCount() - time0) / getTickFrequency();
cout << time0 << endl;
//显示图像
/*double maxVal = 0;
double minVal = 0;
minMaxLoc(mylevelset.m_mPhi, &minVal, &maxVal);
Mat dst = Mat::zeros(src.size(), CV_8U);
mylevelset.m_mPhi.convertTo(dst, CV_8U, 255.0 / (maxVal - minVal), -255.0 * minVal / (maxVal - minVal));
imshow("m_mphi", dst);*/
}
waitKey(0);
return 0;
}
水平集方法：
level_set.h
#pragma once
#include<iostream>
#include<opencv2/opencv.hpp>
#include<string>
class LevelSet
{
public:
LevelSet();
~LevelSet();
//基本参数
int m_iterNum; //迭代次数
float m_lambda1; //全局项系数
float m_nu; //长度约束系数ν
float m_mu; //惩罚项系数μ
float m_timestep; //演化步长δt
float m_epsilon; //规则化参数ε
//过程数据
cv::Mat m_mImage; //源图像
int m_iCol; //图像宽度
int m_iRow; //图像高度
int m_depth; //水平集数据深度
float m_FGValue; //前景值
float m_BKValue; //背景值
//初始化水平集
void initializePhi(cv::Mat img, //输入图像
int iterNum, //迭代次数
cv::Rect boxPhi);//前景区域
void EVolution(cv::String &img_path); //演化
cv::Mat m_mPhi; //水平集：φ
protected:
cv::Mat m_mDirac; //狄拉克处理后水平集：δ（φ）
cv::Mat m_mHeaviside; //海氏函数处理后水平集：Н（φ）
cv::Mat m_mCurv; //水平集曲率κ=div(▽φ/|▽φ|)
cv::Mat m_mK; //惩罚项卷积核
cv::Mat m_mPenalize; //惩罚项中的▽<sup>2</sup>φ
void Dirac(); //狄拉克函数
void Heaviside(); //海氏函数
void Curvature(); //曲率
void BinaryFit(); //计算前景与背景值
};
level_set.cpp
#include"levelset.h"
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.15;//不用调试
m_epsilon = 1.0;
}
LevelSet::~LevelSet()
{
}
void LevelSet::initializePhi(Mat img, int iterNum, Rect boxPhi)
{
//boxPhi是前景区域
m_iterNum = iterNum;
img.copyTo(m_mImage);
//cvtColor(img, m_mImage, COLOR_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, -5.65, 1,
0.5, 1, 0.5);
int c = 15;
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;//在函数外部小于0
}
else
{
if (i>boxPhi.y && i<(boxPhi.y + boxPhi.height) && j>boxPhi.x && j<(boxPhi.x + boxPhi.width))
{
m_mPhi.at<float>(i, j) = c;//在函数内部大于0
}
else
{
m_mPhi.at<float>(i, j) = 0;//在函数上等于0
}
}
}
}
}
//受保护的函数，通过public函数调用
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]);
m_mDirac.at<float>(i, j) = k1 / (k2 + pow(m_mPhi.at<float>(i, j), 2));
}
}
}
//受保护的函数，通过public函数调用
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));
}
}
}
//受保护的函数，通过public函数调用
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;
}
//受保护的函数，通过public函数调用
void LevelSet::BinaryFit()
{
//先计算海氏函数
Heaviside();//得到m_mHeaviside矩阵
//计算前景与背景灰度均值
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;
sumFH += fFHeaviside;
}
}
m_FGValue = sumFG / (sumH + 1e-10); //前景灰度均值
m_BKValue = sumBK / (sumFH + 1e-10); //背景灰度均值
}
//演化函数，调用以上函数
void LevelSet::EVolution(cv::String &img_path)
{
float fCurv;
float fDirac;
float fPenalize;
float fImgValue;
float areaTerm;
float lengthTerm;
float penalizeTerm;
for (int n = 0; n < m_iterNum; n++)
{
Dirac();//得到m_mdirac矩阵
Curvature();//得到m_mcurv矩阵
BinaryFit();//得到m_mHeaviside矩阵，m_FGValue,m_BKValue
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];
lengthTerm = m_nu * fDirac * fCurv; //长度约束
penalizeTerm = m_mu * (fPenalize - fCurv); //惩罚项
//penalizeTerm = m_mu * fDirac;
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);
}
}
//cout << lengthTerm + penalizeTerm + areaTerm << endl;
//显示每一次演化的结果
Mat showIMG = m_mImage.clone();
//cvtColor(m_mImage, showIMG, COLOR_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), 1);
imshow("showIMG", showIMG);
//轮廓中心纵坐标(y)
//轮廓的长轴短轴
//排除重叠血管（lood=0,lood=1）
//相交面积：
waitKey(1);
if (n == 99)
{
//保存文件地址
string save_img_file = "C:\\Users\\tudejiang\\Desktop\\超声图像VOC格式标准数据集201911\\VOC2007\\save1\\";
//保存图片全名
string savedfilename;
int num = img_path.size();
savedfilename = save_img_file + img_path.substr(num - 10);
imwrite(savedfilename,showIMG);
}
}
}
xml解读方法：
官网下载xml解读包
扫一扫
308
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
