Halcon、OpenCV、C++ 實現最小二乘法擬合直線
阿新 • • 發佈：2020-12-28
最小二乘法擬合直線
概念：最小二乘法多項式直線擬合，根據給定的點，求出它的函式y=f(x)，當然求得準確的函式是不太可能的，但是我們能求出它的近似曲線y=φ(x)
原理
假設有點 , I = 1,2,3,……n,求近似曲線y=φ(x)，並且使得y=φ(x)與y=f(x)的平方偏差和最小，偏差
其中我們要找到一組最好的a b ，“最好的”就是要使選出的a b能使得所有的誤差達到最小化。
在此要注意以下，y=ax+b 這裡面的未知量是什麼，很自然的說法是x y,實際上並不是，我們不用去解這個x和y ，因為x和y已經是給定的值了，當我們在找這條直線的時候，我們實際上並不關心x的值有多好，我們要的就是a 和b這兩個變數，它們可以描述x和y之間的關係，我們就是在試圖找出那條最適合的直線所對應的a和b。
可以看到最小二乘法對各個變數求偏導，使得偏導值為0，即可得到最小值，因為e是關於a b的函式，導數為0的點必定是最小值，進入正題
分別對 a b求偏導可以得到：
Halcon最小二乘法擬合直線
1 首先隨機生成一組資料 2 3 Mx:=[100:10:500] 4 5 tuple_length(Mx,len) 6 7 tuple_gen_const(len,5,r) 8 9 Ma:=2 10 11 Mb:=40 12 13 tuple_rand(len , noise) 14 15 My:= Ma *Mx + Mb*noise 16 17 gen_circle(ContCircle, My, Mx, r)
1 接下來用矩陣進行最小二乘求解
2
3 * y = ax + b
4
5 * y1 = ax1 + b
6
7 * y2 = ax2 + b
8
9 * ... .......
10
11 * yn = ax + b
1 create_matrix(len,1,My,y)
2
3 create_matrix(len,2,1,x)
4
5 set_value_matrix(x, [0:len-1], gen_tuple_const(len, 0),Mx)
6
7
8 * XT 代表X的轉置 inv(*)代表*的逆
9
10 * x beta = y
OpenCV最小二乘法擬合直線
1 #include<iostream>
2 #include<opencv2\opencv.hpp>
3 using namespace std;
4 using namespace cv;
5
6
7 int main()
8 {
9 vector<Point>points;
10 //(27 39) (8 5) (8 9) (16 22) (44 71) (35 44) (43 57) (19 24) (27 39) (37 52)
11
12 points.push_back(Point(27, 39));
13 points.push_back(Point(8, 5));
14 points.push_back(Point(8, 9));
15 points.push_back(Point(16, 22));
16 points.push_back(Point(44, 71));
17 points.push_back(Point(35, 44));
18 points.push_back(Point(43, 57));
19 points.push_back(Point(19, 24));
20 points.push_back(Point(27, 39));
21 points.push_back(Point(37, 52));
22 Mat src = Mat::zeros(400, 400, CV_8UC3);
23
24 for (int i = 0;i < points.size();i++)
25 {
26 //在原圖上畫出點
27 circle(src, points[i], 3, Scalar(0, 0, 255), 1, 8);
28 }
29 //構建A矩陣
30 int N = 2;
31 Mat A = Mat::zeros(N, N, CV_64FC1);
32
33 for (int row = 0;row < A.rows;row++)
34 {
35 for (int col = 0; col < A.cols;col++)
36 {
37 for (int k = 0;k < points.size();k++)
38 {
39 A.at<double>(row, col) = A.at<double>(row, col) + pow(points[k].x, row + col);
40 }
41 }
42 }
43 //構建B矩陣
44 Mat B = Mat::zeros(N, 1, CV_64FC1);
45 for (int row = 0;row < B.rows;row++)
46 {
47
48 for (int k = 0;k < points.size();k++)
49 {
50 B.at<double>(row, 0) = B.at<double>(row, 0) + pow(points[k].x, row)*points[k].y;
51 }
52 }
53 //A*X=B
54 Mat X;
55 //cout << A << endl << B << endl;
56 solve(A, B, X,DECOMP_LU);
57 cout << X << endl;
58 vector<Point>lines;
59 for (int x = 0;x < src.size().width;x++)
60 { // y = b + ax;
61 double y = X.at<double>(0, 0) + X.at<double>(1, 0)*x;
62 printf("(%d,%lf)\n", x, y);
63 lines.push_back(Point(x, y));
64 }
65 polylines(src, lines, false, Scalar(255, 0, 0), 1, 8);
66 imshow("src", src);
67
68 //imshow("src", A);
69 waitKey(0);
70 return 0;
71 }
C++最小二乘法擬合直線
1 #include
2 #include
3 #include
4
5 using namespace std;
6
7 int main(int argc, char *argv[])
8 {
9 int num = 0;
10
11 cout << " Input how many numbers you want to calculate:";
12 cin >> num;
13
14 valarray data_x(num);
15 valarray data_y(num);
16
17 while( num )
18 {
19 cout << "Input the "<< num <<" of x:";
20 cin >> data_x[num-1];
21 cout << "Input the "<< num <<" of y:";
22 cin >> data_y[num-1];
23 num--;
24 }
25
26 double A =0.0;
27 double B =0.0;
28 double C =0.0;
29 double D =0.0;
30
31 A = (data_x*data_x).sum();
32 B = data_x.sum();
33 C = (data_x*data_y).sum();
34 D = data_y.sum();
35
36 double k,b,tmp =0;
37 if(tmp=(A*data_x.size()-B*B))
38 {
39 k = (C*data_x.size()-B*D)/tmp;
40 b = (A*D-C*B)/tmp;
41 }
42
43 else
44 {
45 k=1;
46 b=0;
47 }
48
49 cout <<"k="< cout <<"b="<</p>
50
51 return 0;
52 }