多离散点的圆拟合

69 阅读 0 评论 46 点赞

我是靠谱客的博主冷静麦片，最近开发中收集的这篇文章主要介绍多离散点的圆拟合，觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

最近项目涉及到多个圆盘的旋转和运动。这个时候绕不开圆盘圆心、半径的求解。

简单的来说，三点必能确定一个唯一的圆。圆的标准公式是 $（x-x_{0}）^{2}$ $(x-x_0)^2+(y-y_0)^2=r^2$ 。这个公式在求解的时候会比较麻烦，一般会用它的展开公式，就是一般式 $x^2+y^2+Ax+By+C=0$ 。圆心为 $(frac{-A}{2},frac{-B}{2})$ ,圆的半径就是 $sqrt{(frac{-A}{2})^2+(frac{-B}{2})^2-C}$ 。这个通过简单的解方程就可以得到答案了。

下面就是三点求圆的代码：

LONG CalCenter(stLPOSN lTblPosn[] )
{
	stDPOSN	SqPt1, SqPt2, SqPt3 ;
	double	SqPtOneTwo, SqPtOneThree ;
	double	numerator_x, denominator_x, numerator_y, denominator_y ;

	SqPt1.x = pow(lTblPosn[0].x,2) ;
	SqPt1.y = pow(lTblPosn[0].y,2) ;

	SqPt2.x = pow(lTblPosn[1].x,2) ;
	SqPt2.y = pow(lTblPosn[1].y,2) ;

	SqPt3.x = pow(lTblPosn[2].x,2) ;
	SqPt3.y = pow(lTblPosn[2].y,2) ;

	SqPtOneTwo = SqPt1.x + SqPt1.y - SqPt2.x - SqPt2.y ;

	SqPtOneThree = SqPt1.x + SqPt1.y - SqPt3.x - SqPt3.y ;

	numerator_x=2*(lTblPosn[0].y-lTblPosn[1].y)*SqPtOneThree-2*(lTblPosn[0].y-lTblPosn[2].y)*SqPtOneTwo ;

	denominator_x=4*(lTblPosn[0].y-lTblPosn[1].y)*(lTblPosn[0].x-lTblPosn[2].x )-4*lTblPosn[0].y-lTblPosn[2].y)*(lTblPosn[0].x-lTblPosn[1].x) ;
	
	long x = (LONG) ( numerator_x / denominator_x );


	numerator_y =  SqPtOneTwo-2*m_lRingCenter.x*(lTblPosn[0].x- lTblPosn[1].x ) ;
	
	denominator_y =  2 * ( lTblPosn[0].y - lTblPosn[1].y ) ;

	long y = (LONG) ( numerator_y / denominator_y );
	return OK;
}

但实际上，三点算法对于外部输入要求比较高，在很多工程应用中是无法适应的。这个时候就需要找其他算法替代了。一种常用的替代方案就是基于最小二乘法的拟合算法了。最小二乘法是用最简的方法求得一些绝对不可知的真值，而令误差平方之和为最小来寻找一组数据的最佳匹配函数的计算方法，最小二乘法通常用于曲线拟合 (least squares fitting)，拟合算法的推倒过程就省略了，直接看代码吧

int reform_data	(double ax[],
				 double ay[],
				 int	n,
				 double *mx,
				 double *my,
				 double *suu,
				 double *suv,
				 double *svv,
				 double *suuu,
				 double *suvv,
				 double *svvv,
				 double *svuu)
{
	double *au = (double*) malloc(sizeof(double) * n);
	double *av = (double*) malloc(sizeof(double) * n);

	if (au == NULL)
		return -1;		/// not enough memory

	if (av == NULL)
	{
		free(au);
		return -1;
	}
	
	double sx = 0;
	double sy = 0;
	for (int i=0; i<n; i++)
	{
		sx += ax[i];
		sy += ay[i];
	}

	*mx = sx / (double) n;
	*my = sy / (double) n;

	for (int i=0; i<n; i++)
	{
		au[i] = ax[i] - *mx;
		av[i] = ay[i] - *my;
	}


	double uu;
	double uv;
	double vv;

	*suu = 0;
	*suv = 0;
	*svv = 0;
	*suuu = 0;
	*svvv = 0;
	*suvv = 0;
	*svuu = 0;
	for (int i=0; i<n; i++)
	{
		uu = au[i] * au[i];
		vv = av[i] * av[i];
		uv = au[i] * av[i];
		*suu += uu;
		*suv += uv;
		*svv += vv;
		*suuu += au[i] * uu;
		*suvv += uv * av[i];
		*svuu += uv * au[i];
		*svvv += av[i] * vv;
	}

	

	free(au);
	free(av);

	return 0;
}


int fit_circle (double	ax[],
				double	ay[],
				int		n,
				double	*px,
				double	*py,
				double	*pr)
{
	double mx;
	double my;
	double suu;
	double suv;
	double svv;
	double suuu;
	double svuu;
	double suvv;
	double svvv;
	
	int nret = reform_data(ax, ay, n, &mx, &my, &suu, &suv, &svv, &suuu, &suvv, &svvv, &svuu); 
	if (nret != 0)
		return nret;

	double data_a[2][2] = {
		{suu, suv},
		{suv, svv}
	};

	double data_b[] = {
		.5*(suuu+suvv),
		.5*(svvv+svuu)
	};

	double data_x[] = {0, 0};
	
	double det = data_a[0][0]*data_a[1][1] - data_a[1][0] * data_a[0][1];

	data_a[0][0] /= det;
	data_a[0][1] /= det;
	data_a[1][0] /= det;
	data_a[1][1] /= det;

	data_x[0] = data_a[1][1]*data_b[0] - data_a[0][1]*data_b[1];
	data_x[1] = data_a[0][0]*data_b[1] - data_a[1][0]*data_b[0];

	*px = data_x[0]+mx;
	*py = data_x[1]+my;
	*pr = sqrt(data_x[0]*data_x[0]+data_x[1]*data_x[1]+(suu+svv)/n);

	return 0;
}