簡(jiǎn)介：此問題是在做旋轉(zhuǎn)模板匹配的時(shí)候烛愧，選擇最好的匹配結(jié)果時(shí)產(chǎn)生的。查找資料發(fā)現(xiàn)多項(xiàng)式擬合問題可以變成一個(gè)超定方程的求解問題，而opencv中本身有一個(gè)cv::solve()函數(shù)可以求解線性方程組怜姿，因此對(duì)于大多數(shù)用到opencv又要進(jìn)行曲線擬合的地方都可以參考此處的求解過程來解決慎冤。
1. 問題：
原始數(shù)據(jù)是一些離散的散點(diǎn)，下圖是用matplotlib的plot方法畫出來的沧卢，我想得到下圖中最高處附近的近似的曲線方程粪薛，以得到一個(gè)最大值和最大值對(duì)應(yīng)的橫坐標(biāo)。從下圖看搏恤，在最高處附近很像一條拋物線，那就用2次函數(shù)去擬合最高處附近的曲線看看效果
2. 分析
二次函數(shù)的一般形式為,二次函數(shù)由
完全決定湃交，這樣只需要三組
的數(shù)據(jù)我們就可以求出
的表達(dá)式熟空。例如現(xiàn)在我們獲得了三組數(shù)據(jù)，
,寫成方程組的形式就是
求解這個(gè)線性方程組就可以得到我們需要的二次方程搞莺，很容易求出來息罗。
在實(shí)際的情況下我們觀測(cè)獲得的數(shù)據(jù)并不是絕對(duì)準(zhǔn)確的，也就是存在偏差才沧，當(dāng)數(shù)據(jù)足夠多時(shí)就可以去除很大一部分隨機(jī)誤差迈喉，但是當(dāng)方程數(shù)超過未知數(shù)的個(gè)數(shù)時(shí)，怎么求解呢温圆？這就要引入下面超定方程的知識(shí)了挨摸。
3. 超定方程：
設(shè)方程Ax = b.根據(jù)有效的方程個(gè)數(shù)和未知數(shù)的個(gè)數(shù)，可以分為以下3種情況：
rank(A) < n岁歉，也就是說方程個(gè)數(shù)小于未知數(shù)的個(gè)數(shù)得运，約束不夠，方程存在無數(shù)組解锅移，
rank(A) = n 方程個(gè)數(shù)等于未知數(shù)的個(gè)數(shù)熔掺， 方程存在唯一的精確解，解法通常有消元法非剃，LU分解法
rank(A) > n置逻，方程個(gè)數(shù)多于未知數(shù)個(gè)數(shù)，這個(gè)時(shí)候約束過于嚴(yán)格备绽，沒有精確解券坞，這種方程又稱之為超定方程。通常工程應(yīng)用都會(huì)遇到這種情況疯坤，找不到精確解的情況下报慕，選取最優(yōu)解，這個(gè)最優(yōu)解压怠，又稱之為最小二乘解眠冈。
-
超定方程定義：
設(shè)方程組 中，
,
是
維已知向量，
是n維解向量蜗顽，當(dāng)
,即方程個(gè)數(shù)大于自變量個(gè)數(shù)時(shí)布卡，稱此方程組為超定方程組。
- 最小二乘解的定義：
記 ,稱使
最小的解
為方程組
的最小二乘解雇盖。
- 定理：
是
的最小二乘解的充要條件是：
是
的解忿等。
4. 二次曲線擬合：
- 待擬合點(diǎn)：angle為橫坐標(biāo)，matchVal為縱坐標(biāo)崔挖。
angle = { -10., -9., -8., -7., -6., -5., -4., -3., -2., -1., 0.,
1., 2., 3., 4., 5., 6., 7., 8., 9. };
matchVal = { 0.890928, 0.904723, 0.921421, 0.935007, 0.94281 , 0.949828,
0.966265, 0.975411, 0.978693, 0.97662 , 0.974468, 0.967101,
0.957691, 0.958369, 0.949841, 0.932791, 0.9213 , 0.901874,
0.879374, 0.868257 };
- 散點(diǎn)圖：
- 數(shù)學(xué)過程：
設(shè)贸街，則
中
,
,
構(gòu)造，求解此線性方程組就可以得到最小二乘解狸相，也就得到我們需要的二次方程了薛匪。
5. python 實(shí)現(xiàn)：
import cv2 as cv
import numpy as np
import matplotlib.pyplot as plt
#%matplotlib
x = np.array([-10., -9., -8., -7., -6., -5., -4., -3., -2., -1., 0.,
1., 2., 3., 4., 5., 6., 7., 8., 9.], np.float)
y = np.array([0.890928, 0.904723, 0.921421, 0.935007, 0.94281 , 0.949828,
0.966265, 0.975411, 0.978693, 0.97662 , 0.974468, 0.967101,
0.957691, 0.958369, 0.949841, 0.932791, 0.9213 , 0.901874,
0.879374, 0.868257], np.float)
A = np.zeros((len(x), 3)) #構(gòu)造一個(gè)范德蒙德矩陣
A[:,0] = 1
A[:,1] = x
A[:,2] = x**2
c = np.matmul(A.T, A)
d = np.matmul(A.T, y)
_,result = cv.solve(c,d)
y2 = result[0] + result[1]*x + result[2] * (x**2)
plt.grid(True)
plt.xlabel("angle")
plt.ylabel("matchVal")
plt.scatter(x,y)
plt.plot(x,y2)
plt.plot( (-result[1]/result[2]/2,-result[1]/result[2]/2), (0.9,1))
print("擬合方程：f(x) = {} + ({}*x) + ({}*x^2)".format(result[0],result[1],result[2]))
- 輸出結(jié)果：
擬合方程：f(x) = [0.97301156] + ([-0.00231144]*x) + ([-0.00109041]*x^2)
-
擬合結(jié)果：
image.png
6. C++實(shí)現(xiàn)：
#include <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 << "對(duì)稱軸：" << -a1 / a2 / 2 << std::endl;
std::cout << "擬合方程：" << a0 << " + (" << a1 << "x) + (" << a2 << "x^2)";
return 0;
}
int main()
{
double a[] = { -10., -9., -8., -7., -6., -5., -4., -3., -2., -1., 0.,
1., 2., 3., 4., 5., 6., 7., 8., 9. };
double b[] = { 0.890928, 0.904723, 0.921421, 0.935007, 0.94281 , 0.949828,
0.966265, 0.975411, 0.978693, 0.97662 , 0.974468, 0.967101,
0.957691, 0.958369, 0.949841, 0.932791, 0.9213 , 0.901874,
0.879374, 0.868257 };
std::vector<double> x;
std::vector<double> y;
for (int i = 0; i < (sizeof(a) / sizeof(a[0])); i++)
{
x.push_back(a[i]);
y.push_back(b[i]);
}
fitCurve(x, y);
return 0;
}
- 輸出：
A = [1, -10, 100;
1, -9, 81;
1, -8, 64;
1, -7, 49;
1, -6, 36;
1, -5, 25;
1, -4, 16;
1, -3, 9;
1, -2, 4;
1, -1, 1;
1, 0, 0;
1, 1, 1;
1, 2, 4;
1, 3, 9;
1, 4, 16;
1, 5, 25;
1, 6, 36;
1, 7, 49;
1, 8, 64;
1, 9, 81]
b = [0.8909280000000001;
0.9047230000000001;
0.921421;
0.935007;
0.94281;
0.949828;
0.966265;
0.975411;
0.978693;
0.97662;
0.974468;
0.967101;
0.957691;
0.958369;
0.949841;
0.932791;
0.9213;
0.901874;
0.879374;
0.8682569999999999]
result = [0.9730115624060149;
-0.002311438482570065;
-0.001090408407382091]
對(duì)稱軸：-1.0599
擬合方程：0.973012 + (-0.00231144x) + (-0.00109041x^2)
可以和python實(shí)現(xiàn)版對(duì)照著看最后擬合的方程是一樣的！完美脓鹃！
參考：
超定方程理論參考：https://blog.csdn.net/u014652390/article/details/52789591
