OpenCV拟合多项式
在存在多个 (x, y) 点的情况下，要拟合一个多项式 $y = k_0 + k_1x + k_2x^2 + … + k_nx^n$，其中 $n$ 为多项式的阶数。
问题就是求解 $k_0, k_1, k_2, …, k_n$ 的值。
分析
例如对于一个二次多项式 $y = k_0 + k_1x + k_2x^2$，我们可以将其转化为矩阵形式：
\[\begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix} = \begin{bmatrix} 1 & x_0 & x_0^2 \\ 1 & x_1 & x_1^2 \\ 1 & x_2 & x_2^2 \\ \vdots & \vdots & \vdots \\ 1 & x_n & x_n^2 \end{bmatrix} \begin{bmatrix} k_0 \\ k_1 \\ k_2 \end{bmatrix}\]其中， $n$ 是点的个数， $x_i$ 和 $y_i$ 分别是第 $i$ 个点的横坐标和纵坐标， $k_0, k_1, k_2$ 是待求的系数。
使用矩阵来简化表达，记做 $Y = XK$。
要求解这个方程，需要分成三类情况来考虑：
- $\text{rank}(X) < 未知数的个数$。此时方程存在无穷多解。
- $\text{rank}(X) = 未知数的个数$。此时方程存在唯一解。
- $\text{rank}(X) > 未知数的个数$。此时方程没有精确解，只能求出一个近似解。这个最优解被称为最小二乘解。
最小二乘解1.2
对于第三种情况的求解，需要先引入一个概念–超定方程。其定义为：
设方程组$Ax=b$ 中，$A = (a_{ij})_{m \times n}$ ,b是m维已知向量，x是n维解向量，当m> n,即方程个数大于自变量个数时，称此方程组为超定方程组
最小二乘解的定义：
对于超定方程组$Ax=b$，如果存在一个解$x_0$，使得$|Ax_0-b|_2$最小，则称$x_0$为方程组的最小二乘解
求解定理：
设$A = (a_{ij})_{m \times n}$，$b$是$m$维已知向量，$x$是$n$维解向量，$m > n$，则方程组$Ax=b$的最小二乘解为$x_0 = (A^TA)^{-1}A^Tb$
代码实现
C++ 版本1
#include <iostream>
#include <opencv2/opencv.hpp>
int fitCurve(std::vector<double> x, std::vector<double> y)
{
//columns is 3, rows is x.size()
cv::Mat A = cv::Mat::zeros(cv::Size(3, x.size()), CV_64FC1);
for (int i = 0; i < x.size(); i++)
{
A.at<double>(i, 0) = 1;
A.at<double>(i, 1) = x[i];
A.at<double>(i, 2) = x[i] * x[i];
}
cv::Mat b = cv::Mat::zeros(cv::Size(1, y.size()), CV_64FC1);
for (int i = 0; i < y.size(); i++)
{
b.at<double>(i, 0) = y[i];
}
cv::Mat c;
c = A.t() * A;
cv::Mat d;
d = A.t() * b;
cv::Mat result = cv::Mat::zeros(cv::Size(1, 3), CV_64FC1);
cv::solve(c, d, result);
std::cout << "A = " << A << std::endl;
std::cout << "b = " << b << std::endl;
std::cout << "result = " << result << std::endl;
double a0 = result.at<double>(0, 0);
double a1 = result.at<double>(1, 0);
double a2 = result.at<double>(2, 0);
std::cout << "对称轴：" << -a1 / a2 / 2 << std::endl;
std::cout << "拟合方程：" << a0 << " + (" << a1 << "x) + (" << a2 << "x^2)";
return 0;
}
// o(n) 时间复杂度求解 x 对应的 y 值
double polyval(const std::vector<double> k, double x)
{
// vector k is {k0, k1, k2, ..., kn}
double y = 0.;
int degree = k.size() - 1;
for (int i = degree; i >= 0; i--)
{
y = y * x + k[i];
}
return y;
}