CloudCompare——点云二次曲面拟合

本文由CSDN点云侠原创，原文链接。爬虫网站自重。博客长期更新，本文最新更新时间为：2023年8月6日。

1.二次曲面拟合

  二次曲面方程如下：
$$z(x,y)=a{x}^2+bxy+c{y}^2+dx+ey+f\text{\hspace{1em}(1)}$$
采用最小二乘法拟合局部二次曲面，建立目标函数为：
$$z=\sum _{i=0}^k[z(x_i，y_i)-z_i{]}^2\text{\hspace{1em}(2)}$$
联立式( 1) 和式( 2) ，对各项系数求偏导，并令偏导数等于零，得到线性方程组，见式( 3) ，求解线性方程组，得到局部二次拟合曲面方程。
$$\left[\begin{array}{cccccc}{\sum }_{i=0}^k{x}_i^4& {\sum }_{i=0}^k{x}_i^3{y}_i& {\sum }_{i=0}^k{x}_i^2{y}_i& {\sum }_{i=0}^k{x}_i^3& {\sum }_{i=0}^k{x}_i^2{y}_i& {\sum }_{i=0}^k{x}_i^2\\ {\sum }_{i=0}^k{x}_i^3{y}_i& {\sum }_{i=0}^k{x}_i^2{y}_i^2& {\sum }_{i=0}^k{x}_i{y}_i^3& {\sum }_{i=0}^k{x}_i^2{y}_i& {\sum }_{i=0}^k{x}_i{y}_i^2& {\sum }_{i=0}^k{x}_i{y}_i\\ {\sum }_{i=0}^k{x}_i^2{y}_i^2& {\sum }_{i=0}^k{x}_i{y}_i^3& {\sum }_{i=0}^k{y}_i^4& {\sum }_{i=0}^k{x}_i{y}_i^2& {\sum }_{i=0}^k{y}_i^3& {\sum }_{i=0}^k{y}_i^2\\ {\sum }_{i=0}^k{x}_i^3& {\sum }_{i=0}^k{x}_i^2{y}_i& {\sum }_{i=0}^k{x}_i{y}_i^2& {\sum }_{i=0}^k{x}_i^2& {\sum }_{i=0}^k{x}_i{y}_i& {\sum }_{i=0}^k{x}_i\\ {\sum }_{i=0}^k{x}_i^2{y}_i& {\sum }_{i=0}^k{x}_i{y}_i^2& {\sum }_{i=0}^k{y}_i^3& {\sum }_{i=0}^k{x}_i{y}_i& {\sum }_{i=0}^k{y}_i^2& {\sum }_{i=0}^k{y}_i\\ {\sum }_{i=0}^k{x}_i^2& {\sum }_{i=0}^k{x}_i{y}_i& {\sum }_{i=0}^k{y}_i^2& {\sum }_{i=0}^k{x}_i& {\sum }_{i=0}^k{y}_i& {\sum }_{i=0}^{k}1\end{array}\right]\left[\begin{array}{c}a\\ b\\ c\\ d\\ e\\ f\end{array}\right]=\left[\begin{array}{c}{\sum }_{i=0}^k{z}_i{x}_i^2\\ {\sum }_{i=0}^k{z}_i{x}_i{y}_i\\ {\sum }_{i=0}^k{z}_i{y}_i^2\\ {\sum }_{i=0}^k{z}_i{x}_i\\ {\sum }_{i=0}^k{z}_i{y}_i\\ {\sum }_{i=0}^k{z}_i\end{array}\right]\text{\hspace{1em}(3)}$$

2.参考文献

[1] 方程喜,隋立春,朱海雄.用于公路勘测设计的LiDAR点云抽稀算法[J].测绘通报,2017(10):58-61+88.

3.操作流程

  输出二次曲面$a+bx+cy+d{x}^2+exy+f{y}^2=0$的表达式和拟合误差RMS。经实测，CloudCompare的拟合结果偏差较大。

  此外，还生成了拟合曲面的mesh模型，可以用来计算点云到模型的距离。

4.完整操作

5.算法源码

ccQuadric* ccQuadric::Fit(CCLib::GenericIndexedCloudPersist *cloud, double* rms/*=0*/)
{
	//number of points
	unsigned count = cloud->size();
	if (count < CC_LOCAL_MODEL_MIN_SIZE[QUADRIC])
	{
		ccLog::Warning(QString("[ccQuadric::fitTo] Not enough points in input cloud to fit a quadric! (%1 at the very least are required)").arg(CC_LOCAL_MODEL_MIN_SIZE[QUADRIC]));
		return nullptr;
	}

	//project the points on a 2D plane
	CCVector3 G;
	CCVector3 X;
	CCVector3 Y;
	CCVector3 N;
	{
		CCLib::Neighbourhood Yk(cloud);
		
		//plane equation
		const PointCoordinateType* theLSPlane = Yk.getLSPlane();
		if (!theLSPlane)
		{
			ccLog::Warning("[ccQuadric::Fit] Not enough points to fit a quadric!");
			return nullptr;
		}

		assert(Yk.getGravityCenter());
		G = *Yk.getGravityCenter();

		//local base
		N = CCVector3(theLSPlane);
		assert(Yk.getLSPlaneX() && Yk.getLSPlaneY());
		X = *Yk.getLSPlaneX(); //main direction
		Y = *Yk.getLSPlaneY(); //secondary direction
	}

	//project the points in a temporary cloud
	ccPointCloud tempCloud("temporary");
	if (!tempCloud.reserve(count))
	{
		ccLog::Warning("[ccQuadric::Fit] Not enough memory!");
		return nullptr;
	}

	cloud->placeIteratorAtBeginning();
	for (unsigned k=0; k<count; ++k)
	{
		//projection into local 2D plane ref.
		CCVector3 P = *(cloud->getNextPoint()) - G;

		tempCloud.addPoint(CCVector3(P.dot(X),P.dot(Y),P.dot(N)));
	}

	CCLib::Neighbourhood Zk(&tempCloud);
	{
		//set exact values for gravity center and plane equation
		//(just to be sure and to avoid re-computing them)
		Zk.setGravityCenter(CCVector3(0,0,0));
		PointCoordinateType perfectEq[4] = { 0, 0, 1, 0 };
		Zk.setLSPlane(	perfectEq,
						CCVector3(1,0,0),
						CCVector3(0,1,0),
						CCVector3(0,0,1));
	}

	Tuple3ub dims;
	const PointCoordinateType* eq = Zk.getQuadric(&dims);
	if (!eq)
	{
		ccLog::Warning("[ccQuadric::Fit] Failed to fit a quadric!");
		return nullptr;
	}

	//we recenter the quadric object
	ccGLMatrix glMat(X,Y,N,G);

	ccBBox bb = tempCloud.getOwnBB();
	CCVector2 minXY(bb.minCorner().x,bb.minCorner().y);
	CCVector2 maxXY(bb.maxCorner().x,bb.maxCorner().y);

	ccQuadric* quadric = new ccQuadric(minXY, maxXY, eq, &dims, &glMat);

	quadric->setMetaData(QString("Equation"),QVariant(quadric->getEquationString()));

	//compute rms if necessary
	if (rms)
	{
		const unsigned char dX = dims.x;
		const unsigned char dY = dims.y;
		//const unsigned char dZ = dims.z;

		*rms = 0;

		for (unsigned k=0; k<count; ++k)
		{
			//projection into local 2D plane ref.
			const CCVector3* P = tempCloud.getPoint(k);

			PointCoordinateType z = eq[0] + eq[1]*P->u[dX] + eq[2]*P->u[dY] + eq[3]*P->u[dX]*P->u[dX] + eq[4]*P->u[dX]*P->u[dY] + eq[5]*P->u[dY]*P->u[dY];
			*rms += (z-P->z)*(z-P->z);
		}

		if (count)
		{
			*rms = sqrt(*rms / count);
			quadric->setMetaData(QString("RMS"),QVariant(*rms));
		}
	}
	
	return quadric;
}

6.相关代码

[1] PCL 最小二乘拟合二次曲面
[2] matlab 最小二乘拟合二次曲面
[3] Open3D 点云最小二乘拟合二次曲面