OpenCV4开发入门教程099:曲线拟合

索引地址:系列索引

数学上的曲线拟合使用的是最小二乘法。

原理

最小二乘法多项式曲线拟合,根据给定的m个点,并不要求这条曲线精确地经过这些点,而是曲线y=f(x)的近似曲线y=φ(x)y= \varphi(x)

给定数据点pi(xi,yi)p_i(x_i,y_i),其中i=1,2,…,m。求近似曲线y=φ(x)y= \varphi(x)。并且使得近似曲线与y=f(x)的偏差最小。近似曲线在点pip_i处的偏差δi=φ(xi)yi=1,2,...,m\delta i= \varphi(x_i)-y,i=1,2,...,m

常见的曲线拟合方法

使偏差绝对值之和最小

minφi=1mδi=i=1mφ(xi)yi\min\limits_{\varphi} \sum\limits_{i=1}^{m}|\delta_i|=\sum\limits_{i=1}^{m}|\varphi(x_i)-y_i|

使偏差绝对值最大的最小

minφmaxiδi=φ(xi)yi\min\limits_{\varphi} \max\limits_{i}|\delta_{i}|=|\varphi(x_i)-y_i|

使偏差平方和最小

minφi=1mδi2=i=1m(φ(xi)yi)2\min\limits_{\varphi} \sum\limits_{i=1}^{m}\delta_{i}^{2}=\sum\limits_{i=1}^{m}(\varphi(x_i)-y_i)^2

按偏差平方和最小的原则选取拟合曲线,并且采取二项式方程为拟合曲线的方法,称为最小二乘法。

推到过程

1.设拟合多项式为

y=a0+a1x+...+akxky=a_0+a_1 x + ... + a_k x^k

2.各点到这条曲线的距离之和,即偏差平方和如下

R2=i=1n[yi(a0+a1xi+...+akxik)]2R^2=\sum\limits_{i=1}^{n}[y_i-(a_0+a_1 x_i + ... + a_k x_i^k)]^2

3.为了求得符合条件的a值,对等式右边求ai偏导数,因而我们得到了

2i=1n[y(a0+a1x+...+akxk)]x=02i=1n[y(a0+a1x+...+akxk)]=0...2i=1n[y(a0+a1x+...+akxk)]xk=0-2\sum\limits_{i=1}^{n}[y-(a_0 + a_1 x + ... + a_k x^k)]x=0 \\ -2\sum\limits_{i=1}^{n}[y-(a_0+a_1 x + ... + a_k x^k)]=0 \\ ... \\ -2\sum\limits_{i=1}^{n}[y-(a_0 + a_1 x + ... + a_k x^k)]x^k=0

4.将等式左边进行一下化简,然后应该可以得到下面的等式

a0n+a1i=1nxi+...+aki=1nxika0i=1nxi+a1i=1nxi2+...+aki=1nxik+1...a0i=1nxik+a1i=1nxik+1+...+aki=1nxi2ka_0 n + a_1 \sum\limits_{i=1}^{n} x_i + ... + a_k \sum\limits_{i=1}^{n}x_i^k \\ a_0\sum\limits_{i=1}^{n}x_i + a_1 \sum\limits_{i=1}^{n}x_i^2 + ... + a_k \sum\limits_{i=1}^{n}x_i^{k+1} \\ ... \\ a_0 \sum\limits_{i=1}^{n} x_i^k + a_1 \sum\limits_{i=1}^{n} x_i^{k+1} + ... + a_k \sum\limits_{i=1}^{n}x_i^{2k}

5.把这些等式表示成矩阵的形式,就可以得到下面的矩阵

[ni=1nxi...i=1nxiki=1nxii=1nxi2...i=1nxik+1.........i=1nxiki=1nxik+1...i=1nxi2k][a0a1...ak]=[i=1nyii=1nxiyi...i=1nxikyi]\left[ \begin{array}{l} n & \sum\limits_{i=1}^{n}x_i & ... & \sum\limits_{i=1}^{n}x_i^k \\ \sum\limits_{i=1}^{n}x_i & \sum\limits_{i=1}^{n}x_i^2 & ... & \sum\limits_{i=1}^{n}x_i^{k+1} \\ . & . & . \\ . & . & . \\ . & . & . \\ \sum\limits_{i=1}^{n}x_i^k & \sum\limits_{i=1}^{n}x_i^{k+1} & ... & \sum\limits_{i=1}^{n}x_i^{2k} \end{array} \right] \left[\begin{array}{l} a_0 \\ a_1 \\ . \\ . \\ . \\ a_k \end{array}\right]= \left[\begin{array}{l} \sum\limits_{i=1}^{n}y_i \\ \sum\limits_{i=1}^{n}x_iy_i \\ . \\ . \\ . \\ \sum\limits_{i=1}^{n}x_i^k y_i \end{array}\right]

6.将这个范德蒙得矩阵化简后可得到

[1x1...x1k1x2...x2k............1xn...xnk][a0a1...ak]=[y1y2...yn]\left[\begin{array}{l} 1 & x_1 & ... & x_1^k \\ 1 & x_2 & ... & x_2^k \\ . & . & . & . \\ . & . & . & . \\ . & . & . & . \\ 1 & x_n & ... & x_n^k \end{array}\right] \left[\begin{array}{l} a_0 \\ a_1 \\ . \\ . \\ . \\ a_k \end{array}\right]= \left[\begin{array}{l} y_1 \\ y_2 \\ . \\ . \\ . \\ y_n \end{array}\right]

7.也就是说XA=YX*A=Y,那么A=(XX)1XYA = (X'*X)-1*X'*Y,便得到了系数矩阵A,同时,我们也就得到了拟合曲线。

在opencv中,有一个专门用于求解线性方程的函数,即cv::solve(),具体调用形式如下:

1
2
3
4
5
6
int cv::solve(
cv::InputArray X, // 左边矩阵X, nxn
cv::InputArray Y, // 右边矩阵Y,nx1
cv::OutputArray A, // 结果,系数矩阵A,nx1
int method = cv::DECOMP_LU // 估算方法
);

我们只需要按照上述原理,构造出矩阵X和Y,即可调用该函数,计算出多项式的系数矩阵A。

opencv中支持的估算方法如下所示:

  • DECOMP_LU 选择最佳枢轴元素的高斯消去法。
  • DECOMP_SVD 奇异值分解(SVD)方法;系统可能被过度定义和/或矩阵src1可能是奇异的
  • DECOMP_EIG 特征值分解;矩阵src1必须是对称的
  • DECOMP_CHOLESKY Cholesky LLTLL^T分解;矩阵src1必须是对称且正定义的
  • DECOMP_QR QR分解;系统可能被过度定义和/或矩阵src1可能是奇异的
  • DECOMP_NORMAL 虽然所有之前的标志都是互斥的,但该标志可以与任何之前的标志一起使用;这意味着普通方程src1Tsrc1dst=src1Tsrc2src1^T⋅src1⋅dst=src1^Tsrc2被求解,而不是原始系统中的src1⋅dst=src2

像销售额、出货量、访问量等等都可以使用这个方法进行预测。使用已有数据拟合出数据方程,然后得出将来的预测值。

测试代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
bool polynomial_curve_fit(std::vector<cv::Point>& key_point, int n, cv::Mat& A)
{
//Number of key points
int N = key_point.size();

//构造矩阵X
cv::Mat X = cv::Mat::zeros(n + 1, n + 1, CV_64FC1);
for (int i = 0; i < n + 1; i++)
{
for (int j = 0; j < n + 1; j++)
{
for (int k = 0; k < N; k++)
{
X.at<double>(i, j) = X.at<double>(i, j) +
std::pow(key_point[k].x, i + j);
}
}
}

//构造矩阵Y
cv::Mat Y = cv::Mat::zeros(n + 1, 1, CV_64FC1);
for (int i = 0; i < n + 1; i++)
{
for (int k = 0; k < N; k++)
{
Y.at<double>(i, 0) = Y.at<double>(i, 0) +
std::pow(key_point[k].x, i) * key_point[k].y;
}
}

A = cv::Mat::zeros(n + 1, 1, CV_64FC1);
//求解矩阵A
cv::solve(X, Y, A, cv::DECOMP_LU); //X*A = Y
return true;
}

int main()
{
//创建用于绘制的深蓝色背景图像
cv::Mat image = cv::Mat::zeros(480, 640, CV_8UC3);
image.setTo(cv::Scalar(100, 0, 0));

//输入拟合点
std::vector<cv::Point> points;
points.push_back(cv::Point(100., 58.));
points.push_back(cv::Point(150., 70.));
points.push_back(cv::Point(200., 90.));
points.push_back(cv::Point(252., 140.));
points.push_back(cv::Point(300., 220.));
points.push_back(cv::Point(350., 400.));

//将拟合点绘制到空白图上
for (int i = 0; i < points.size(); i++)
{
cv::circle(image, points[i], 5, cv::Scalar(0, 0, 255), 2, 8, 0);
}

//绘制折线
cv::polylines(image, points, false, cv::Scalar(0, 255, 0), 1, 8, 0);

cv::Mat A;

polynomial_curve_fit(points, 3, A);
std::cout << "A = " << A << std::endl;

std::vector<cv::Point> points_fitted;

for (int x = 0; x < 400; x++)
{
double y = A.at<double>(0, 0) + A.at<double>(1, 0) * x +
A.at<double>(2, 0)*std::pow(x, 2) + A.at<double>(3, 0)*std::pow(x, 3);

points_fitted.push_back(cv::Point(x, y));
}
cv::polylines(image, points_fitted, false, cv::Scalar(0, 255, 255), 1, 8, 0);

cv::imshow("image", image);

cv::waitKey(0);
return 0;
}

输出为

1
2
3
4
A = [-80.59143874478204;
2.591187861030039;
-0.01564864659563801;
3.472543637058225e-05]

效果为

result

OpenCV4开发入门教程099:曲线拟合
https://feater.top/opencv/opencv-outline-fit
作者
JackeyLea
发布于
2020年10月6日
许可协议