OpenCV3.3最小二乘法直线拟合

// LineFit.cpp : 此文件包含 "main" 函数。程序执行将在此处开始并结束。 // #include "pch.h" #include <iostream> #include <opencv2/opencv.hpp> #include <opencv2/highgui/highgui.hpp> #include <opencv2/imgproc/imgproc.hpp> using namespace std; using namespace cv; //直线拟合 static void LinearFitting(vector<double> x,vector<double> y, double &slope,double &intercept,double &r_square) { int size_x = x.size(); int size_y = y.size(); double xmean = 0.0; double ymean = 0.0; for (size_t i = 0; i < size_x; i++) { xmean += x.at(i); ymean += y.at(i); } xmean /= size_x; ymean /= size_y; double sumx2 = 0.0; double sumy2 = 0.0; double sumxy = 0.0; for (size_t i = 0; i < size_x; i++) { sumx2 += (x.at(i) - xmean)*(x.at(i) - xmean); sumy2 += (y.at(i) - ymean)*(y.at(i) - ymean); sumxy += (y.at(i) - ymean)*(x.at(i) - xmean); } slope = sumxy / sumx2; intercept = ymean - slope * xmean; r_square = sumxy * sumxy / (sumx2*sumy2); } 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 (size_t 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; } // 运行程序: Ctrl + F5 或调试 >“开始执行(不调试)”菜单 // 调试程序: F5 或调试 >“开始调试”菜单 // 入门提示: // 1. 使用解决方案资源管理器窗口添加/管理文件 // 2. 使用团队资源管理器窗口连接到源代码管理 // 3. 使用输出窗口查看生成输出和其他消息 // 4. 使用错误列表窗口查看错误 // 5. 转到“项目”>“添加新项”以创建新的代码文件，或转到“项目”>“添加现有项”以将现有代码文件添加到项目 // 6. 将来，若要再次打开此项目，请转到“文件”>“打开”>“项目”并选择 .sln 文件