Halcon、OpenCV、C++ 实现最小二乘法拟合直线

最小二乘法拟合直线

概念:最小二乘法多项式直线拟合,根据给定的点,求出它的函数y=f(x),当然求得准确的函数是不太可能的,但是我们能求出它的近似曲线y=φ(x)

 


原理

假设有点  , I = 1,2,3,……n,求近似曲线y=φ(x),并且使得y=φ(x)与y=f(x)的平方偏差和最小,偏差

image.png

image.png

其中我们要找到一组最好的a b ,“最好的”就是要使选出的a b能使得所有的误差达到最小化。


在此要注意以下,y=ax+b 这里面的未知量是什么,很自然的说法是x y,实际上并不是,我们不用去解这个x和y ,因为x和y已经是给定的值了,当我们在找这条直线的时候,我们实际上并不关心x的值有多好,我们要的就是a 和b这两个变量,它们可以描述x和y之间的关系,我们就是在试图找出那条最适合的直线所对应的a和b。

image.png


image.png

可以看到最小二乘法对各个变量求偏导,使得偏导值为0,即可得到最小值,因为e是关于a b的函数,导数为0的点必定是最小值,进入正题

 分别对 a b求偏导可以得到:

image.png

image.png

Halcon最小二乘法拟合直线

*首先随机生成一组数据

Mx:=[100:10:500]
tuple_length(Mx,len)
tuple_gen_const(len,5,r)
Ma:=2
Mb:=40
tuple_rand(len , noise)
My:= Ma *Mx + Mb*noise
gen_circle(ContCircle, My, Mx, r)

image.png

接下来用矩阵进行最小二乘求解

* y = ax + b

* y1 = ax1 + b

* y2 = ax2 + b

* ... .......

* yn = ax + b

image.png

create_matrix(len,1,My,y)
create_matrix(len,2,1,x)
set_value_matrix(x, [0:len-1], gen_tuple_const(len, 0),Mx)


* XT 代表X的转置 inv(*)代表*的逆
* x beta = y
* xT x beta = xT y
* beta = inv( xT x) xT y
mult_matrix(x,x,'ATB',xtx)
mult_matrix(x,y,'ATB',xty)
invert_matrix(xtx,'general', 0, invxtx)
mult_matrix(invxtx,xty,'AB', beta)
get_full_matrix(beta, Values)
Newy:=Values[0] * [10,800] + Values[1]
gen_contour_polygon_xld(Contour, Newy, [10,800])

image.png

OpenCV最小二乘法拟合直线

关于solve函数 ,可以在这个链接查看用法:solve函数使用

#include<iostream>
#include<opencv2\opencv.hpp>
using namespace std;
using namespace cv;
 
 
int main()
{
	vector<Point>points;
	//(27 39) (8 5) (8 9) (16 22) (44 71) (35 44) (43 57) (19 24) (27 39) (37 52)
 
	points.push_back(Point(27, 39));
	points.push_back(Point(8, 5));
	points.push_back(Point(8, 9));
	points.push_back(Point(16, 22));
	points.push_back(Point(44, 71));
	points.push_back(Point(35, 44));
	points.push_back(Point(43, 57));
	points.push_back(Point(19, 24));
	points.push_back(Point(27, 39));
	points.push_back(Point(37, 52));
	Mat src = Mat::zeros(400, 400, CV_8UC3);
 
	for (int i = 0;i < points.size();i++)
	{
		//在原图上画出点
		circle(src, points[i], 3, Scalar(0, 0, 255), 1, 8);
	}
	//构建A矩阵 
	int N = 2;
	Mat A = Mat::zeros(N, N, CV_64FC1);
 
	for (int row = 0;row < A.rows;row++)
	{
		for (int col = 0; col < A.cols;col++)
		{
			for (int k = 0;k < points.size();k++)
			{
				A.at<double>(row, col) = A.at<double>(row, col) + pow(points[k].x, row + col);
			}
		}
	}
	//构建B矩阵
	Mat B = Mat::zeros(N, 1, CV_64FC1);
	for (int row = 0;row < B.rows;row++)
	{
 
		for (int k = 0;k < points.size();k++)
		{
			B.at<double>(row, 0) = B.at<double>(row, 0) + pow(points[k].x, row)*points[k].y;
		}
	}
	//A*X=B
	Mat X;
	//cout << A << endl << B << endl;
	solve(A, B, X,DECOMP_LU);
	cout << X << endl;
	vector<Point>lines;
	for (int x = 0;x < src.size().width;x++)
	{				// y = b + ax;
		double y = X.at<double>(0, 0) + X.at<double>(1, 0)*x;
		printf("(%d,%lf)\n", x, y);
		lines.push_back(Point(x, y));
	}
	polylines(src, lines, false, Scalar(255, 0, 0), 1, 8);
	imshow("src", src);
	
	//imshow("src", A);
	waitKey(0);
	return 0;
}

image.png

C++最小二乘法拟合直线

#include 
#include 
#include
 
using namespace std;
 
int main(int argc, char *argv[])
{
int num = 0;
 
cout << " Input how many numbers you want to calculate:";
cin >> num;
 
valarray data_x(num);
valarray data_y(num);
 
while( num )
{
cout << "Input the "<< num <<" of x:";
cin >> data_x[num-1];
cout << "Input the "<< num <<" of y:";
cin >> data_y[num-1];
num--;
}
 
double A =0.0;
double B =0.0;
double C =0.0;
double D =0.0;
 
A = (data_x*data_x).sum();
B = data_x.sum();
C = (data_x*data_y).sum();
D = data_y.sum();
 
double k,b,tmp =0;
if(tmp=(A*data_x.size()-B*B))
{
k = (C*data_x.size()-B*D)/tmp;
b = (A*D-C*B)/tmp;
}
 
else
{
k=1;
b=0;
}
 
cout <<"k="< cout <<"b="<</p>
 
return 0;
}


————————————————

版权声明:本文为CSDN博主「c1learning」的原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接及本声明。

原文链接:https://blog.csdn.net/c1learning/article/details/101447218


本文出自勇哥的网站《少有人走的路》wwww.skcircle.com,转载请注明出处!讨论可扫码加群:

发表评论:

◎欢迎参与讨论,请在这里发表您的看法、交流您的观点。

会员中心
搜索
«    2024年4月    »
1234567
891011121314
15161718192021
22232425262728
2930
网站分类
标签列表
最新留言
    热门文章 | 热评文章 | 随机文章
文章归档
友情链接
  • 订阅本站的 RSS 2.0 新闻聚合
  • 扫描加本站机器视觉QQ群,验证答案为:halcon勇哥的机器视觉
  • 点击查阅微信群二维码
  • 扫描加勇哥的非标自动化群,验证答案:C#/C++/VB勇哥的非标自动化群
  • 扫描加站长微信:站长微信:abc496103864
  • 扫描加站长QQ:
  • 扫描赞赏本站:
  • 留言板:

Powered By Z-BlogPHP 1.7.2

Copyright Your skcircle.com Rights Reserved.

鄂ICP备18008319号


站长QQ:496103864 微信:abc496103864