diff --git a/sourceCode/dataFitting.cpp b/sourceCode/dataFitting.cpp new file mode 100644 index 0000000..33dab6a --- /dev/null +++ b/sourceCode/dataFitting.cpp @@ -0,0 +1,385 @@ +#include "SG_baseDataType.h" +#include "SG_baseAlgo_Export.h" +#include +#ifdef __WIN32 +#include +#endif // __WIN32 + +#include +#include +#include + +void lineFitting(std::vector< SVzNL3DPoint>& inliers, double* _k, double* _b) +{ + //最小二乘拟合直线参数 + double xx_sum = 0; + double x_sum = 0; + double y_sum = 0; + double xy_sum = 0; + int num = 0; + for (int i = 0; i < inliers.size(); i++) + { + x_sum += inliers[i].x; //x的累加和 + y_sum += inliers[i].y; //y的累加和 + xx_sum += inliers[i].x * inliers[i].x; //x的平方累加和 + xy_sum += inliers[i].x * inliers[i].y; //x,y的累加和 + num++; + } + *_k = (num * xy_sum - x_sum * y_sum) / (num * xx_sum - x_sum * x_sum); //根据公式求解k + *_b = (-x_sum * xy_sum + xx_sum * y_sum) / (num * xx_sum - x_sum * x_sum);//根据公式求解b +} + +//拟合成通用直线方程ax+by+c=0,包括垂直 +void lineFitting_abc(std::vector< SVzNL3DPoint>& inliers, double* _a, double* _b, double* _c) +{ + //判断是否为垂直 + int dataSize = (int)inliers.size(); + if (dataSize < 2) + return; + + double deltaX = abs(inliers[0].x - inliers[dataSize - 1].x); + double deltaY = abs(inliers[0].y - inliers[dataSize - 1].y); + std::vector< SVzNL3DPoint> fittingData; + if (deltaX < deltaY) + { + //x=ky+b 拟合 + for (int i = 0; i < dataSize; i++) + { + SVzNL3DPoint a_fitPt; + a_fitPt.x = inliers[i].y; + a_fitPt.y = inliers[i].x; + a_fitPt.z = inliers[i].z; + fittingData.push_back(a_fitPt); + } + double k = 0, b = 0; + lineFitting(fittingData, &k, &b); + //ax+by+c + *_a = 1.0; + *_b = -k; + *_c = -b; + } + else + { + //y = kx+b拟合 + double k = 0, b = 0; + lineFitting(inliers, &k, &b); + //ax+by+c + *_a = k; + *_b = -1; + *_c = b; + } + return; +} + +//圆最小二乘拟合 +double fitCircleByLeastSquare( + const std::vector& pointArray, + SVzNL3DPoint& center, + double& radius) +{ + int N = pointArray.size(); + if (N < 3) { + return std::numeric_limits::max(); + } + + double sumX = 0.0; + double sumY = 0.0; + double sumX2 = 0.0; + double sumY2 = 0.0; + double sumX3 = 0.0; + double sumY3 = 0.0; + double sumXY = 0.0; + double sumXY2 = 0.0; + double sumX2Y = 0.0; + + for (int pId = 0; pId < N; ++pId) { + sumX += pointArray[pId].x; + sumY += pointArray[pId].y; + + double x2 = pointArray[pId].x * pointArray[pId].x; + double y2 = pointArray[pId].y * pointArray[pId].y; + sumX2 += x2; + sumY2 += y2; + + sumX3 += x2 * pointArray[pId].x; + sumY3 += y2 * pointArray[pId].y; + sumXY += pointArray[pId].x * pointArray[pId].y; + sumXY2 += pointArray[pId].x * y2; + sumX2Y += x2 * pointArray[pId].y; + } + + double C, D, E, G, H; + double a, b, c; + + C = N * sumX2 - sumX * sumX; + D = N * sumXY - sumX * sumY; + E = N * sumX3 + N * sumXY2 - (sumX2 + sumY2) * sumX; + G = N * sumY2 - sumY * sumY; + H = N * sumX2Y + N * sumY3 - (sumX2 + sumY2) * sumY; + + a = (H * D - E * G) / (C * G - D * D); + b = (H * C - E * D) / (D * D - G * C); + c = -(a * sumX + b * sumY + sumX2 + sumY2) / N; + + center.x = -a / 2.0; + center.y = -b / 2.0; + radius = sqrt(a * a + b * b - 4 * c) / 2.0; + + double err = 0.0; + double e; + double r2 = radius * radius; + for (int pId = 0; pId < N; ++pId) { + e = pow(pointArray[pId].x - center.x, 2) + pow(pointArray[pId].y - center.y, 2) - r2; + if (e > err) { + err = e; + } + } + return err; +} + +#if 0 +bool leastSquareParabolaFit(const std::vector& points, + double& a, double& b, double& c, + double& mse, double& max_err) +{ + // 校验点集数量(至少3个点才能拟合抛物线) + int n = points.size(); + if (n < 3) { + return false; + } + + // 初始化各项求和参数 + double sum_x = 0.0, sum_x2 = 0.0, sum_x3 = 0.0, sum_x4 = 0.0; + double sum_y = 0.0, sum_xy = 0.0, sum_x2y = 0.0; + + // 计算各项求和 + for (const auto& p : points) { + double x = p.x; + double y = p.y; + double x2 = x * x; + double x3 = x2 * x; + double x4 = x3 * x; + + sum_x += x; + sum_x2 += x2; + sum_x3 += x3; + sum_x4 += x4; + sum_y += y; + sum_xy += x * y; + sum_x2y += x2 * y; + } + + // 构建线性方程组 M * [a,b,c]^T = N + // M矩阵:3x3 + double M[3][3] = { + {sum_x4, sum_x3, sum_x2}, + {sum_x3, sum_x2, sum_x}, + {sum_x2, sum_x, (double)n} + }; + + // N向量:3x1 + double N[3] = { sum_x2y, sum_xy, sum_y }; + + // 高斯消元法求解线性方程组(3元一次方程组) + // 步骤1:将M转化为上三角矩阵 + for (int i = 0; i < 3; i++) { + // 选主元(避免除数为0) + int pivot = i; + for (int j = i; j < 3; j++) { + if (fabs(M[j][i]) > fabs(M[pivot][i])) { + pivot = j; + } + } + // 交换行 + std::swap(M[i], M[pivot]); + std::swap(N[i], N[pivot]); + + // 归一化主元行 + double div = M[i][i]; + if (fabs(div) < 1e-10) { + return false; + } + for (int j = i; j < 3; j++) { + M[i][j] /= div; + } + N[i] /= div; + + // 消去下方行 + for (int j = i + 1; j < 3; j++) { + double factor = M[j][i]; + for (int k = i; k < 3; k++) { + M[j][k] -= factor * M[i][k]; + } + N[j] -= factor * N[i]; + } + } + + // 步骤2:回代求解 + c = N[2]; + b = N[1] - M[1][2] * c; + a = N[0] - M[0][1] * b - M[0][2] * c; + + // 计算拟合误差 + mse = 0.0; + max_err = 0.0; + for (const auto& p : points) { + double y_fit = a * p.x * p.x + b * p.x + c; + double err = y_fit - p.y; + double err_abs = fabs(err); + mse += err * err; + if (err_abs > max_err) { + max_err = err_abs; + } + } + mse /= n; // 均方误差 + + return true; +} +#endif +//抛物线最小二乘拟合 y=ax^2 + bx + c +bool leastSquareParabolaFitEigen( + const std::vector& points, + double& a, double& b, double& c, + double& mse, double& max_err) +{ + int n = points.size(); + if (n < 3) { + return false; + } + + // 构建系数矩阵A和目标向量B + Eigen::MatrixXd A(n, 3); + Eigen::VectorXd B(n); + for (int i = 0; i < n; i++) { + double x = points[i].x; + A(i, 0) = x * x; + A(i, 1) = x; + A(i, 2) = 1.0; + B(i) = points[i].y; + } + + // 最小二乘求解:Ax = B,直接调用Eigen的求解器 + Eigen::Vector3d coeffs = A.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(B); + a = coeffs(0); + b = coeffs(1); + c = coeffs(2); + + // 计算误差 + mse = 0.0; + max_err = 0.0; + for (const auto& p : points) { + double y_fit = a * p.x * p.x + b * p.x + c; + double err = y_fit - p.y; + double err_abs = fabs(err); + mse += err * err; + if (err_abs > max_err) { + max_err = err_abs; + } + } + mse /= n; + return true; +} + +//计算面参数: z = Ax + By + C +//res: [0]=A, [1]= B, [2]=-1.0, [3]=C, +void vzCaculateLaserPlane(std::vector Points3ds, std::vector& res) +{ + //最小二乘法拟合平面 + //获取cv::Mat的坐标系以纵向为x轴,横向为y轴,而cvPoint等则相反 + //系数矩阵 + cv::Mat A = cv::Mat::zeros(3, 3, CV_64FC1); + // + cv::Mat B = cv::Mat::zeros(3, 1, CV_64FC1); + //结果矩阵 + cv::Mat X = cv::Mat::zeros(3, 1, CV_64FC1); + double x2 = 0, xiyi = 0, xi = 0, yi = 0, zixi = 0, ziyi = 0, zi = 0, y2 = 0; + for (int i = 0; i < Points3ds.size(); i++) + { + x2 += (double)Points3ds[i].x * (double)Points3ds[i].x; + y2 += (double)Points3ds[i].y * (double)Points3ds[i].y; + xiyi += (double)Points3ds[i].x * (double)Points3ds[i].y; + xi += (double)Points3ds[i].x; + yi += (double)Points3ds[i].y; + zixi += (double)Points3ds[i].z * (double)Points3ds[i].x; + ziyi += (double)Points3ds[i].z * (double)Points3ds[i].y; + zi += (double)Points3ds[i].z; + } + A.at(0, 0) = x2; + A.at(1, 0) = xiyi; + A.at(2, 0) = xi; + A.at(0, 1) = xiyi; + A.at(1, 1) = y2; + A.at(2, 1) = yi; + A.at(0, 2) = xi; + A.at(1, 2) = yi; + A.at(2, 2) = (double)((int)Points3ds.size()); + B.at(0, 0) = zixi; + B.at(1, 0) = ziyi; + B.at(2, 0) = zi; + //计算平面系数 + X = A.inv() * B; + //A + res.push_back(X.at(0, 0)); + //B + res.push_back(X.at(1, 0)); + //Z的系数为-1 + res.push_back(-1.0); + //C + res.push_back(X.at(2, 0)); + return; +} + +/** + * @brief 空间直线最小二乘拟合 + * @param points 输入的三维点集(至少2个点,否则无意义) + * @param center 输出:拟合直线的质心(基准点P0) + * @param direction 输出:拟合直线的方向向量v(单位化) + * @return 拟合是否成功(点集有效返回true) + */ +bool fitLine3DLeastSquares(const std::vector& points, SVzNL3DPoint& center, SVzNL3DPoint& direction) +{ + // 检查点集有效性 + if (points.size() < 2) { + std::cerr << "Error: 点集数量必须大于等于2!" << std::endl; + return false; + } + + int n = points.size(); + Eigen::MatrixXd A(n, 3); // 点集矩阵:每行一个点的(x,y,z) + + // 1. 计算质心(center) + double cx = 0.0, cy = 0.0, cz = 0.0; + for (const auto& p : points) { + cx += p.x; + cy += p.y; + cz += p.z; + A.row(points.size() - n) << p.x, p.y, p.z; // 填充点集矩阵 + n--; + } + cx /= points.size(); + cy /= points.size(); + cz /= points.size(); + center = { cx, cy, cz }; + + // 2. 构造去中心化的协方差矩阵(3x3) + // 关键修复:使用RowVector3d(行向量)做rowwise减法,匹配维度 + Eigen::RowVector3d centroid_row(cx, cy, cz); + Eigen::MatrixXd centered = A.rowwise() - centroid_row; // 维度匹配,无报错 + + // 协方差矩阵计算(n-1为无偏估计,工程中也可直接用n) + Eigen::Matrix3d cov = centered.transpose() * centered; // / (points.size() - 1); + // 3. 特征值分解:求协方差矩阵的特征值和特征向量 + Eigen::SelfAdjointEigenSolver eigensolver(cov); + if (eigensolver.info() != Eigen::Success) { + std::cerr << "Error: 特征值分解失败!" << std::endl; + return false; + } + + // 最大特征值对应的特征向量即为方向向量(Eigen默认按特征值升序排列,取最后一个) + Eigen::Vector3d dir = eigensolver.eigenvectors().col(2); + // 单位化方向向量(可选,但工程中通常标准化) + dir.normalize(); + + direction = { dir(0), dir(1), dir(2) }; + return true; +}