SVD拟合直线、平面

最新推荐文章于 2025-03-03 07:30:00 发布

原创最新推荐文章于 2025-03-03 07:30:00 发布 · 330 阅读

1 ·

CC 4.0 BY-SA版权

文章标签：

#c++ #几何学 #算法 #机器学习

SVD的原理、应用参考博客：奇异值分解（SVD）方法求解最小二乘问题的原理 - 一抹烟霞 - 博客园

SVD拟合直线：SVD的第一列就是直线的方向，直线位置点通过中心化测量点集得到

SVD拟合平面：SVD第一列向量与第二列向量的叉乘就是平面的法向，平面的基点通过法向上的平移得到

具体步骤如下（以下代码使用了GSL数学库和opencascade模型库，且仅使用直线点集的中心点作为位置点）：

	// 中心化测量点，并计算中心点作为拟合直线的位置。
		gp_Pnt location;
		const vector<gp_Pnt> points = Centralize(measuredPoints, location);

		// 创建矩阵。
		gsl_matrix* gslA = gsl_matrix_alloc(points.size(), 3);
		unique_ptr<gsl_matrix, function<void(gsl_matrix*)>> gslA_ForAutoDelete(gslA, GslMatrixDeleter);

		for (int i = 0; i < points.size(); i++) {
			gsl_matrix_set(gslA, i, 0, points[i].X());
			gsl_matrix_set(gslA, i, 1, points[i].Y());
			gsl_matrix_set(gslA, i, 2, points[i].Z());
		}

		// SVD分解。
		gsl_matrix* gslV = gsl_matrix_alloc(3, 3);
		unique_ptr<gsl_matrix, function<void(gsl_matrix*)>> gslV_ForAutoDelete(gslV, GslMatrixDeleter);

		gsl_vector* gslD = gsl_vector_alloc(3);
		unique_ptr<gsl_vector, function<void(gsl_vector*)>> gslD_ForAutoDelete(gslD, GslVectorDeleter);

		gsl_vector* work = gsl_vector_alloc(3);
		unique_ptr<gsl_vector, function<void(gsl_vector*)>> work_ForAutoDelete(work, GslVectorDeleter);

		int gsl_errno = gsl_linalg_SV_decomp(gslA, gslV, gslD, work);
		if (gsl_errno) {
			return nullptr;
		}

		// 第一个右奇异向量作为直线方向。
		gsl_vector* gslDirection = gsl_vector_alloc(3);
		unique_ptr<gsl_vector, function<void(gsl_vector*)>> gslDirection_ForAutoDelete(gslDirection, GslVectorDeleter);
		gsl_matrix_get_col(gslDirection, gslV, 0);

		gp_XYZ gpDirectinXYZ;
		GslFacility::Copy(gslDirection, gpDirectinXYZ);

		// 得到拟合直线。
		return make_shared<gp_Ax1>(location, gpDirectinXYZ);