利用opencv进行高斯函数拟合
opencv
OpenCV: 开源计算机视觉库
项目地址：https://gitcode.com/gh_mirrors/opencv31/opencv
免费下载资源
·
摘要
在实际应用中，很多数据都符合高斯分布。比如正态分布就是高斯分布，γ光子和电子发生碰撞反应后产生的两个电子所具有的能量值也是符合高斯分布的。而所谓的高斯函数拟合就是利用大量符合高斯分布的数据点进行拟合，从而得到具体的高斯函数。本文只研究一维高斯函数的拟合。
一维高斯函数基础
高斯函数，Gaussian Function， 也简称为Gaussian，一维形式如下：
对于任意的实数a,b,c，是以著名数学家Carl Friedrich Gauss的名字命名的。高斯的一维图是特征对称“bell curve”形状，a是曲线尖峰的高度，b是尖峰中心的坐标，c称为标准方差，表征的是bell钟状的宽度。注意：参数a，b，c的位置无所谓。
高斯函数曲线如下图所示：
一维高斯函数拟合原理
所谓的高斯函数拟合就是利用大量二维数据点(x1,y1),(x2,y2),(x3,y3), … 来求解出参数a，b，c。显然直接求解难度很大，所以我们转换为多项式拟合。如下所述：
设有一组实验数据(x_{i},y_{i})(i = 1,2,3,…N),可用高斯函数描述：
式（3.1）中待估参数a,c和b，分别代表的物理意义为高斯曲线的峰高、峰位置和半宽度信息。将式(3.1)两边取自然对数，化为：
令：
则式（3.2）化为二次多项式拟合函数：
考虑全部数据和测量误差，并以矩阵形式表示如下：
简记为：
在不考虑总量程误差E的影响情况下，根据最小二乘原理，可求得拟合常数b0，b1，b2构成的矩阵B的广义最小二乘解为：
进而根据式（3.3），可求得待估参数a，b，c：
将所得的a，b和c带入式（3.1），就能够得到由实验数据(x_{i},y_{i})(i = 1,2,3,…N)拟合的高斯函数。
最小二乘法
最小二乘法（又称最小平方法）是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。最小二乘法还可用于曲线拟合。其他一些优化问题也可通过最小化能量或最大化熵用最小二乘法来表达。
1801年，意大利天文学家朱赛普·皮亚齐发现了第一颗小行星谷神星。经过40天的跟踪观测后，由于谷神星运行至太阳背后，使得皮亚齐失去了谷神星的位置。随后全世界的科学家利用皮亚齐的观测数据开始寻找谷神星，但是根据大多数人计算的结果来寻找谷神星都没有结果。时年24岁的高斯也计算了谷神星的轨道。奥地利天文学家海因里希·奥尔伯斯根据高斯计算出来的轨道重新发现了谷神星。
高斯使用的最小二乘法的方法发表于1809年他的著作《天体运动论》中。法国科学家勒让德于1806年独立发明“最小二乘法”，但因不为世人所知而默默无闻。勒让德曾与高斯为谁最早创立最小二乘法原理发生争执。1829年，高斯提供了最小二乘法的优化效果强于其他方法的证明，因此被称为高斯-马尔可夫定理。
opencv进行高斯拟合
根据前面的知识我们可以知道，只要求解出矩阵B就可以求解出高斯待估参数a，b，c。而在opencv中，有一个专门用于求解线性方程的函数，即cv::solve()，具体调用形式如下：
我们只需要按照上述原理，构造出矩阵X和Z，即可调用该函数，从而计算出多项式的系数矩阵B。
opencv中支持的估算方法如下图所示：
实现代码
首先我们需要构造矩阵X和Z，假设数据点保存在 vectorcv::Point& cvFitPointVec 这个容器中，矩阵X如下所示：
所以矩阵X的构造代码如下所示：
int N = cvFitPointVec.size();
// n就是矩阵X的cols(列数)，即 n+1 = 3，所以n=2
// N 是二维数据点的个数
cv::Mat X(N, n + 1, CV_32FC1, cv::Scalar(0));
for (int i = 0; i < N; ++i)
{
float* ptr = X.ptr<float>(i);
for (int j = 0; j < n + 1; ++j)
{
ptr[j] = pow(cvFitPointVec[i].x, j);
}
}
从 式3.3 可知 矩阵Z的元素值为 lnyi (i = 0,1,2,…) ，为了方便我们将矩阵Z改名为矩阵Y，所以矩阵Y构造的伪代码如下所示：
int N = cvFitPointVec.size();
cv::Mat Y(N, 1, CV_32FC1, cv::Scalar(0));
for (int i = 0; i < N; ++i)
{
Y.at<float>(i, 0) = log(cvFitPointVec[i].y);
}
实现代码如下所示：
bool gaussiantCurveFit(const vector<cv::Point>& cvFitPointVec, int n, cv::Mat& A)
{
int N = cvFitPointVec.size();
cv::Mat Y(N, 1, CV_32FC1, cv::Scalar(0));
cv::Mat X(N, n + 1, CV_32FC1, cv::Scalar(0));
for (int i = 0; i < N; ++i)
{
float* ptr = X.ptr<float>(i);
for (int j = 0; j < n + 1; ++j)
{
ptr[j] = pow(cvFitPointVec[i].x, j);
}
Y.at<float>(i, 0) = log(cvFitPointVec[i].y);
}
//求解矩阵A：X*A = Y
solve(X, Y, A, cv::DECOMP_QR);
return true;
}
测试代码
为了便于测试，这里我们提供了符合高斯分布的二维数据点，保存在csv表格中，所以我们需要线读取csv表格，将数据存储到 vectorcv::PointcvFitPointVec 这个容器中，然后就可以调用 gaussiantCurveFit 函数得到线性矩阵B，从而计算出高斯参数a，b，c。最后我们在同一个窗口中用不同颜色绘制出csv表格中二维数据点的连线以及拟合后的高斯函数曲线。
利用csv的表格功能，我们可以初略查看出数据点的分布，如下图所示：
//原始数据点
vector < cv::Point> g_cvRawPointVec;
//拟合数据点
vector < cv::Point> g_fitPointVec;
const cv::Size kImgLabelSize = Size(500, 150);
void readCsv()
{
ifstream inFile("./data.csv",ios::in);
if (!inFile)
{
cout << "打开文件失败！" << endl;
exit(1);
}
g_cvRawPointVec.reserve(3000);
string line;
cv::Point point;
while (getline(inFile, line))
{
string field;
istringstream sin(line);
getline(sin, field, ',');
//将刚刚读取的字符串转换成int
int x = atoi(field.c_str());
getline(sin, field, ',');
//将刚刚读取的字符串转换成int
int y = atoi(field.c_str());
getline(sin, field);
g_cvRawPointVec.push_back(cv::Point(x, y));
}
inFile.close();
}
bool gaussiantCurveFit(const vector<cv::Point>& cvRawPointVec, int n, cv::Mat& A)
{
int N = cvRawPointVec.size();
cv::Mat Y(N, 1, CV_32FC1, cv::Scalar(0));
cv::Mat X(N, n + 1, CV_32FC1, cv::Scalar(0));
for (int i = 0; i < N; ++i)
{
float* ptr = X.ptr<float>(i);
for (int j = 0; j < n + 1; ++j)
{
ptr[j] = pow(cvRawPointVec[i].x, j);
}
Y.at<float>(i, 0) = log(cvRawPointVec[i].y);
}
//求解矩阵A：X*A = Y
solve(X, Y, A, cv::DECOMP_QR);
return true;
}
void getMat(cv::Mat& mat, float fParmA, float fParmB, float fParmC)
{
uint32_t uiXMaxValue = 3000;
uint32_t uiYMaxValue = 200;
double dXAsixRatio = (double)(uiXMaxValue - 1) / (double)kImgLabelSize.width;
double dYAsixRatio;
kImgLabelSize.height > (int)uiYMaxValue ? dYAsixRatio = kImgLabelSize.height / uiYMaxValue : dYAsixRatio = uiYMaxValue / kImgLabelSize.height;
for (size_t i = 0; i < g_cvRawPointVec.size(); ++i)
{
if (kImgLabelSize.height > uiYMaxValue)
{
g_cvRawPointVec[i] = cv::Point(g_cvRawPointVec[i].x / dXAsixRatio, kImgLabelSize.height - g_cvRawPointVec[i].y * dYAsixRatio);
}
else
{
g_cvRawPointVec[i] = cv::Point(g_cvRawPointVec[i].x / dXAsixRatio, kImgLabelSize.height - g_cvRawPointVec[i].y / dYAsixRatio);
}
}
//原始数据点的折线图
polylines(mat, g_cvRawPointVec, false, cv::Scalar(255, 0, 0), 1, 8, 0);
g_fitPointVec.reserve(uiXMaxValue);
double dMaxYValue = -1.0;
for (uint32_t x = 0; x < uiXMaxValue; ++x)
{
double y = 0.0;
if ((int)fParmB != 0)
{
y = fParmA * exp(-pow(x - fParmC, 2) / fParmB);
}
if (y >= dMaxYValue)
{
dMaxYValue = y;
}
g_fitPointVec.push_back(cv::Point(x, y));
}
dXAsixRatio = (uiXMaxValue - 1) / (double)kImgLabelSize.width;
kImgLabelSize.height > uiYMaxValue ? dYAsixRatio = kImgLabelSize.height / uiYMaxValue : dYAsixRatio = uiYMaxValue / kImgLabelSize.height;
for (uint32_t i = 0; i < uiXMaxValue; ++i)
{
if (kImgLabelSize.height > uiYMaxValue)
{
g_fitPointVec[i] = cv::Point(g_fitPointVec[i].x / dXAsixRatio, kImgLabelSize.height - g_fitPointVec[i].y * dYAsixRatio);
}
else
{
g_fitPointVec[i] = cv::Point(g_fitPointVec[i].x / dXAsixRatio, kImgLabelSize.height - g_fitPointVec[i].y / dYAsixRatio);
}
}
polylines(mat, g_fitPointVec, false, cv::Scalar(0, 0, 255), 1, 8, 0);
}
int main()
{
readCsv();
cv::Mat A(3, 1, CV_32FC1, cv::Scalar(0));
gaussiantCurveFit(g_cvRawPointVec, 2, A);
float a0 = A.at<float>(0, 0);
float a1 = A.at<float>(1, 0);
float a2 = A.at<float>(2, 0);
float fParmB = -1 / (a2);
float fParmC = -a1 / (2 * a2);
float fParmA = std::exp(a0 + fParmC * fParmC / fParmB);
cv::Mat mat = cv::Mat::zeros(kImgLabelSize.height, kImgLabelSize.width, CV_8UC3);
mat.setTo(cv::Scalar(255, 255, 255));
getMat(mat,fParmA,fParmB,fParmC);
imshow("mat", mat);
moveWindow("mat", 500, 400);
waitKey(0);
return 0;
}
运行结果：
途中红色曲线是拟合后的高斯曲线，蓝色是原始数据点。从图像上看拟合的效果不是很好，上述代码中我们将csv中所有数据都进行拟合，但是这样效果并不好。为了提高拟合效果，我们并不将所有数据点都进行拟合，而是先找到纵坐标最大的点，然后以这个点为中心，左右各选取100个数据点，这样总共选取了201个数据点进行拟合。理论上说这样的拟合应该是有效果的。
修改后的代码：
//原始数据点
vector<cv::Point> g_cvRawPointVec;
//拟合数据点
vector<cv::Point> g_fitPointVec;
const cv::Size kImgLabelSize = Size(500, 150);
cv::Point g_cvPhotopeakPoint = cv::Point(-1, -1);
void readCsv()
{
ifstream inFile("./data.csv",ios::in);
if (!inFile)
{
cout << "打开文件失败！" << endl;
exit(1);
}
g_cvRawPointVec.reserve(3000);
string line;
cv::Point point;
while (getline(inFile, line))//getline(inFile, line)表示按行读取CSV文件中的数据
{
string field;
//将整行字符串line读入到字符串流sin中
istringstream sin(line);
//将字符串流sin中的字符读入到field字符串中，以逗号为分隔符
getline(sin, field, ',');
//将刚刚读取的字符串转换成int
int x = atoi(field.c_str());
//将字符串流sin中的字符读入到field字符串中，以逗号为分隔符
getline(sin, field, ',');
//将刚刚读取的字符串转换成int
int y = atoi(field.c_str());
getline(sin, field); //将字符串流sin中的字符读入到field字符串中，以逗号为分隔符
g_cvRawPointVec.push_back(cv::Point(x, y));
}
inFile.close();
}
void setFitPointVec()
{
for (int i = 0; i < g_cvRawPointVec.size(); ++i)
{
if (g_cvRawPointVec[i].y >= g_cvPhotopeakPoint.y)
{
g_cvPhotopeakPoint = g_cvRawPointVec[i];
}
}
if (g_cvRawPointVec.size() > 201)
{
auto itr = find(g_cvRawPointVec.begin(), g_cvRawPointVec.end(), g_cvPhotopeakPoint);
int iMinSpace = 0;
itr - g_cvRawPointVec.begin() > g_cvRawPointVec.end() - itr ?
iMinSpace = g_cvRawPointVec.end() - itr :
iMinSpace = itr - g_cvRawPointVec.begin();
if (iMinSpace < 201)
{
g_fitPointVec.assign(g_cvRawPointVec.begin(), g_cvRawPointVec.end());
}
else
{
if (itr + 201 == g_cvRawPointVec.end())
{
g_fitPointVec.assign(itr - 201, itr + 201);
}
else
{
g_fitPointVec.assign(itr - 201, itr + 201 + 1);
}
}
}
else
{
g_fitPointVec.assign(g_cvRawPointVec.begin(), g_cvRawPointVec.end());
}
}
bool gaussiantCurveFit(const vector<cv::Point>& cvRawPointVec, int n, cv::Mat& A)
{
int N = cvRawPointVec.size();
cv::Mat Y(N, 1, CV_32FC1, cv::Scalar(0));
cv::Mat X(N, n + 1, CV_32FC1, cv::Scalar(0));
for (int i = 0; i < N; ++i)
{
float* ptr = X.ptr<float>(i);
for (int j = 0; j < n + 1; ++j)
{
ptr[j] = pow(cvRawPointVec[i].x, j);
}
Y.at<float>(i, 0) = log(cvRawPointVec[i].y);
}
//求解矩阵A：X*A = Y
solve(X, Y, A, cv::DECOMP_QR);
return true;
}
void getMat(cv::Mat& mat, float fParmA, float fParmB, float fParmC)
{
uint32_t uiXMaxValue = 3000;
uint32_t uiYMaxValue = 200;
double dXAsixRatio = (double)(uiXMaxValue - 1) / (double)kImgLabelSize.width;
double dYAsixRatio;
kImgLabelSize.height > (int)uiYMaxValue ? dYAsixRatio = kImgLabelSize.height / uiYMaxValue : dYAsixRatio = uiYMaxValue / kImgLabelSize.height;
for (size_t i = 0; i < g_cvRawPointVec.size(); ++i)
{
if (kImgLabelSize.height > uiYMaxValue)
{
//由于后面没有用到m_cvRawPointVec，所以直接修改m_cvRawPointVec
g_cvRawPointVec[i] = cv::Point(g_cvRawPointVec[i].x / dXAsixRatio, kImgLabelSize.height - g_cvRawPointVec[i].y * dYAsixRatio);
}
else
{
g_cvRawPointVec[i] = cv::Point(g_cvRawPointVec[i].x / dXAsixRatio, kImgLabelSize.height - g_cvRawPointVec[i].y / dYAsixRatio);
}
}
//原始数据点的折线图
polylines(mat, g_cvRawPointVec, false, cv::Scalar(255, 0, 0), 1, 8, 0);
vector<cv::Point> cvFitPoint;
cvFitPoint.reserve(uiXMaxValue);
double dMaxYValue = -1.0;
for (uint32_t x = 0; x < uiXMaxValue; ++x)
{
double y = 0.0;
if ((int)fParmB != 0)
{
y = fParmA * exp(-pow(x - fParmC, 2) / fParmB);
}
if (y >= dMaxYValue)
{
dMaxYValue = y;
}
cvFitPoint.push_back(cv::Point(x, y));
}
dXAsixRatio = (uiXMaxValue - 1) / (double)kImgLabelSize.width;
kImgLabelSize.height > uiYMaxValue ? dYAsixRatio = kImgLabelSize.height / uiYMaxValue : dYAsixRatio = uiYMaxValue / kImgLabelSize.height;
for (uint32_t i = 0; i < uiXMaxValue; ++i)
{
if (kImgLabelSize.height > uiYMaxValue)
{
cvFitPoint[i] = cv::Point(cvFitPoint[i].x / dXAsixRatio, kImgLabelSize.height - cvFitPoint[i].y * dYAsixRatio);
}
else
{
cvFitPoint[i] = cv::Point(cvFitPoint[i].x / dXAsixRatio, kImgLabelSize.height - cvFitPoint[i].y / dYAsixRatio);
}
}
polylines(mat, cvFitPoint, false, cv::Scalar(0, 0, 255), 1, 8, 0);
}
int main()
{
readCsv();
setFitPointVec();
cv::Mat A(3, 1, CV_32FC1, cv::Scalar(0));
gaussiantCurveFit(g_fitPointVec, 2, A);
float a0 = A.at<float>(0, 0);
float a1 = A.at<float>(1, 0);
float a2 = A.at<float>(2, 0);
float fParmB = -1 / (a2);
float fParmC = -a1 / (2 * a2);
float fParmA = std::exp(a0 + fParmC * fParmC / fParmB);
cv::Mat mat = cv::Mat::zeros(kImgLabelSize.height, kImgLabelSize.width, CV_8UC3);
mat.setTo(cv::Scalar(255, 255, 255));
getMat(mat,fParmA,fParmB,fParmC);
imshow("mat", mat);
moveWindow("mat", 500, 400);
waitKey(0);
return 0;
}
最终运行结果：
从上面效果图我们可以看到拟合的效果很好，完全符合高斯分布的趋势。至此利用opencv进行高斯曲线拟合的相关知识全部书写完毕。由于不能上传csv表格，所以不能提供demo中的数据。
GitHub 加速计划 / opencv31 / opencv
162
15
下载
OpenCV: 开源计算机视觉库
最近提交(Master分支：4 个月前 )
796adf5d
Avoid adding value to nullptr 14 小时前
1a6ef7e0
Fixed Android build with Vulkan support. 17 小时前
旨在为数千万中国开发者提供一个无缝且高效的云端环境，以支持学习、使用和贡献开源项目。
更多推荐
4
0- 0
已为社区贡献1条内容
所有评论(0)