C++ 最小二乘法 直线拟合、曲线拟合、平面拟合、高斯拟合

本节介绍如何用Eigen求解线性最小二乘方程组。求解Ax=b的最小二乘问题，等价于求解方程

 使用Eigen的求解的代码如下：

    Eigen::MatrixXd MatX;//样本数据Eigen::MatrixXd MatY;//观测值Eigen::MatrixXd MatLS;//待定系数MatLS = (MatX.transpose() * MatX).inverse() * MatX.transpose() * MatY;

 直线拟合：

直线方程：y = k*x + b

//y = k*x + b
void least_square_line(double input[][2], int number, double&k, double&b){Eigen::MatrixXd MatX;Eigen::MatrixXd MatY;MatX.resize(number, 2);MatY.resize(number, 1);Eigen::MatrixXd MatLS;MatLS.resize(number, 1);for(int i = 0; i < number; ++i){MatX.row(i) = Eigen::Vector2d(input[i][0], 1);MatY(i, 0) = input[i][1];}MatLS = (MatX.transpose()*MatX).inverse()*MatX.transpose()*MatY;//MatLS = MatX.colPivHouseholderQr().solve(MatY); //也可以用k = MatLS(0,0);b = MatLS(1,0);
}

曲线拟合：

 多项式方程：Pn(x)=a(n)x^n+a(n-1)x^(n-1)+…+a(1)x+a(0)

//y = a0 + a1*x + a2*x*x + a3*x*x*x + a4*x*x*x*x + a5*x*x*x*x*x;
void least_square_curve(double input[][2], int number, int jie, std::vector&weight){Eigen::MatrixXd MatX;Eigen::MatrixXd MatY;MatX.resize(number, jie);MatY.resize(number, 1);Eigen::MatrixXd MatLS;MatLS.resize(number, 1);for(int i = 0; i < number; ++i){std::vector pows;for(int j = 0; j < jie; j++){double p =  pow(input[i][0], j);MatX(i, j) = p;}MatY.row(i)[0] = input[i][1];}MatLS = (MatX.transpose()*MatX).inverse()*MatX.transpose()*MatY;for(int i = 0; i < jie; i++)weight.push_back(MatLS(i, 0)) ;
}

3阶多项式

5阶多项式：

7阶多项式：

平面拟合： 

 平面方程：Ax+By+Cz+D=0

//z = A*x + B*y + C
void Least_squares_plane(double input[][3], int number, double& A, double& B, double& C)
{Eigen::MatrixXd MatX;Eigen::MatrixXd MatY;MatX.resize(number, 3);MatY.resize(number, 1);Eigen::MatrixXd MatLS;MatLS.resize(number, 1);for (int i = 0; i < number; ++i) {MatX.row(i) = Eigen::RowVector3d(input[i][0], input[i][1],1);MatY.row(i)[0] = input[i][2];}MatLS = (MatX.transpose() * MatX).inverse() * MatX.transpose() * MatY;A = MatLS(0, 0);B = MatLS(1, 0);C = MatLS(2, 0);}

以该3D模型的每个三角形面片的顶点为观测数据：

拟合的平面如下：

高斯拟合：

 高斯函数

//ln(z) = A + B*x + C*y + D*x*x + E*y*y
//z = exp(A + B*x + C*y + D*x*x + E*y*y)
void Least_squares_gaussian(double input[][3], int number, double& A, double& B, double& C, double& D, double& E)
{Eigen::MatrixXd MatX;Eigen::MatrixXd MatY;MatX.resize(number, 5);MatY.resize(number, 1);Eigen::MatrixXd MatLS;MatLS.resize(number, 1);for (int i = 0; i < number; ++i) {Eigen::Matrix row;row << 1 , input[i][0] , input[i][1] , input[i][0] * input[i][0] , input[i][1] * input[i][1];MatX.row(i) = row;MatY.row(i)[0] = log(input[i][2]);}MatLS = (MatX.transpose() * MatX).inverse() * MatX.transpose() * MatY;A = MatLS(0, 0);B = MatLS(1, 0);C = MatLS(2, 0);D = MatLS(3, 0);E = MatLS(4, 0);
}

同样以该3D模型的顶点为观测数据

拟合高斯面如下：

本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！