图像处理中的椭圆拟合(一)

        图像处理中的椭圆检测用处还是挺多的,找到这里来的同学大多是想用椭圆检测来解决某些实际问题吧,所以我就不做介绍,直奔主题。我研究这块也有一段时间了,也查找了挺多资料,貌似通用的椭圆算法还没有,不排除我孤陋寡闻了。前辈提出的算法适用范围比较有限,这个比较有限是相对直线检测来说的。但直接用Hough变换来找椭圆几乎是不可能的事,在5维空间里做投票,想想都觉得可怕。于是有人想到用随机Hough变换。这是一种很合理的方法,我就是这么做的,不过这种方法有个不足之处,后面会讲到。这里先介绍这方法的流程。

        二次曲线的一般方程为:,其中(x,y)为图像的坐标空间,BCDEF是二次曲线的参数。当满足时二次曲线为椭圆。方程中需要求解的参数有5个,在随机Hough变换过程中要至少采集5个点,得到5个方程以求得这5个参数。若方程有解且满足约束条件,则将解加入参数空间进行累积。思路是比较简单的,下面边贴代码边解释(P.S. 代码仅供参考)。

void FindEllipse(TImage* OrgGrayImg)
{
int ImgHeight = OrgGrayImg->nHeight;
	int ImgWidth = OrgGrayImg->nWidth;
	unsigned char * Img = OrgGrayImg->pImage; // 输入图像确保是二值图像

srand((unsigned)time(NULL));

	int totalPt = 0;// 用于统计样本点的个数

	for (i = 0; i < ImgHeight - 0; i++)
	{
		unsigned char *imgdata = Img + i * ImgWidth;
		for (j = 0; j < ImgWidth - 0; j++)
		{
			if (!imgdata[j])
				totalPt ++;
		}
	}

	if (totalPt < 5)
		return;

	POINT * seq;
	seq = new POINT [totalPt];

	int count = 0;
	for (i = 0; i < ImgHeight; i++)
	{
		unsigned char *data = Img + i * ImgWidth;
		for (j = 0; j < ImgWidth; j++)
		{
			if (!data[j])
			{
				seq[count].x = j;
				seq[count].y = i;
				count ++;
			}
		}
	}

	double Para[5];	// 存放结果(5个参数A,B,C,D,E)的数组
	int Angle_V[360]={0};// 椭圆倾斜角参数空间
	int *Center_XV = new int[ImgWidth];// 椭圆中心点x坐标参数空间
	int *Center_YV = new int[ImgHeight];// 椭圆中心点y坐标参数空间
	int *A_axis_V = new int[max(ImgWidth,ImgHeight)/2];// 椭圆长轴参数空间
	int *B_axis_V = new int[max(ImgWidth,ImgHeight)/2];// 椭圆短轴参数空间
	
	memset(Center_XV,0,sizeof(int)*ImgWidth);
	memset(Center_YV,0,sizeof(int)*ImgHeight);
	memset(A_axis_V,0,sizeof(int)*max(ImgWidth,ImgHeight)/2);
	memset(B_axis_V,0,sizeof(int)*max(ImgWidth,ImgHeight)/2);

	double Theta,X_c,Y_c,A_axis,B_axis;

int loop = 1;// 成功求出参数的迭代次数
	int looptop = loop * 1;// 总的迭代次数(也就是控制计算时间的上限,以免陷入无限循环)
while(loop > 0 && looptop > 0)
{
looptop --;
	int idx;
	for (count = totalPt; count > 0; count--)// 打乱样本点排列的顺序
	{
		POINT ptrtmp;
		idx = rand() % count;
		
		ptrtmp = seq[idx];
		seq[idx] = seq[count-1];
		seq[count-1] = ptrtmp;		
	}

	double PioMatrix[5*5];
	for (i = 0; i < 5; i++)
	{
		PioMatrix[i*5] = seq[i].x * seq[i].x;
		PioMatrix[i*5 + 1] = 2 * seq[i].x * seq[i].y;
		PioMatrix[i*5 + 2] = seq[i].y * seq[i].y;
		PioMatrix[i*5 + 3] = 2 * seq[i].x;
		PioMatrix[i*5 + 4] = 2 * seq[i].y;
	}

	if (GaussJordanInv(PioMatrix,5) == false)// Gauss-Jordan求逆
		continue;
	double sum;
	for (i = 0; i < 5; i++)
	{
		sum = 0;
		for (j = 0; j < 5; j++)
		{
			sum +=  -(PioMatrix[i*5 + j]);
		}
		Para[i] = sum;
	}

	if (pow(Para[1],2) - Para[0] * Para[2] > 0)
			continue;

		if (fabs(Para[0] - Para[2]) < 1e-20)
			Theta = 1.570796326;
		else if (Para[0] > Para[2])
			Theta = 0.5 * (atan(2.0 * Para[1] / (Para[0] - Para[2])) + PI);
		else
			Theta = 0.5 * (atan(2.0 * Para[1] / (Para[0] - Para[2])));

		X_c = (4.0 * Para[1] * Para[4] - 4.0 * Para[2] * Para[3]) / (4.0 * Para[0] * Para[2] - 4.0 * Para[1] * Para[1]);
		Y_c = (4.0 * Para[1] * Para[3] - 4.0 * Para[0] * Para[4]) / (4.0 * Para[0] * Para[2] - 4.0 * Para[1] * Para[1]);
		A_axis = 2 * (Para[0] * pow(X_c,2) + Para[2] * pow(Y_c,2) + 2 * Para[1] * X_c * Y_c - 1) 
						/ (Para[0] + Para[2] - sqrt(pow(Para[0] - Para[2],2) + pow(2.0 * Para[1],2)));
		B_axis = 2 * (Para[0] * pow(X_c,2) + Para[2] * pow(Y_c,2) + 2 * Para[1] * X_c * Y_c - 1) 
						/ (Para[0] + Para[2] + sqrt(pow(Para[0] - Para[2],2) + pow(2.0 * Para[1],2)));
		
		A_axis = sqrt(A_axis);	//长轴
		B_axis = sqrt(B_axis);	//短轴

		int AngleTmp = (int)(Theta * 180 / PI + 360 + 0.5) % 360;
		Angle_V[AngleTmp]++;
		
		if (X_c < 0 || Y_c < 0 || A_axis < 0 || B_axis < 0)
			continue;
		if (X_c >= ImgWidth || Y_c >= ImgHeight || A_axis > max(ImgWidth,ImgHeight)/2 || B_axis > max(ImgWidth,ImgHeight)/2)
			continue;

		if (X_c >= 0 && X_c < ImgWidth)
			Center_XV[(int)X_c]++;
		if (Y_c >= 0 && Y_c < ImgHeight)
			Center_YV[(int)Y_c]++;
		if (A_axis >= 0 && A_axis < max(ImgWidth,ImgHeight)/2)
			A_axis_V[(int)A_axis]++;
		if (B_axis >= 0 && B_axis < max(ImgWidth,ImgHeight)/2)
			B_axis_V[(int)B_axis]++;		
		loop--;
}
	
	int Angle,Ai,Bi,Cx,Cy;
	//	Angle
	int MaxPara = 0;
	for (i = 0; i < 360; i++)
	{
		if (MaxPara < Angle_V[i])
		{
			MaxPara = Angle_V[i];
			Angle = i;
		}
	}
	//	Cy
	MaxPara = 0;
	for (i = 0; i < ImgHeight; i++)
	{
		if (MaxPara < Center_YV[i])
		{
			MaxPara = Center_YV[i];
			Cy = i;
		}
	}
	//	Cx
	MaxPara = 0;
	for (i = 0; i < ImgWidth; i++)
	{
		if (MaxPara < Center_XV[i])
		{
			MaxPara = Center_XV[i];
			Cx = i;
		}
	}
	//	Ai
	MaxPara = 0;
	for (i = 0; i < max(ImgWidth,ImgHeight)/2; i++)
	{
		if (MaxPara < A_axis_V[i])
		{
			MaxPara = A_axis_V[i];
			Ai = i;
		}
	}
	//	Bi
	MaxPara = 0;
	for (i = 0; i < max(ImgWidth,ImgHeight)/2; i++)
	{
		if (MaxPara < B_axis_V[i])
		{
			MaxPara = B_axis_V[i];
			Bi = i;
		}
	}

	delete[] Center_XV;
	delete[] Center_YV;
	delete[] A_axis_V;
	delete[] B_axis_V;


	double sma = SinMem[Angle];
	double cma = CosMem[Angle];
	for (int n = 0; n < 360; n++)
	{
		i = (Bi) * CosMem[360 - n];
		j = (Ai) * SinMem[360 - n];
		
		int x,y;
		x = (j * cma - i * sma) + Cx;
		y = (i * cma + j * sma) + Cy;
		
		Mask[y * ImgWidth + x] = 0;
	}
	delete[] Mask;
	delete[] seq;
}

测试结果:

原图:

拟合结果(虚线为拟合的椭圆):

前面说到这种方法有缺陷,请看下面的情形:

原图:

拟合结果:

       当样本点只集中在椭圆的一边时,随机5点的hough变换总会拟合错误,实际应用中往往会发生这样的情况。这是因为公式错了吗?于是我单独提取出五点做测试,即只做一次迭代。测试结果如下图所示。图中实线为实际椭圆,我是用画图工具拖出来的“完美椭圆”,用橡皮擦擦掉一大半部分,最后再做一点旋转。打交叉的是取样的5点,虚线是用这5点代入公式求得的拟合椭圆。可见求得的椭圆穿过了5个点,表明不是求解错了,可就是跟实际的有很大差别。唯一的解释是取样点不在我们想要的椭圆上,也就是说即使是用画图工具拖出来看似完美的椭圆并不完美,这是因为样本点的坐标是整型,精度很低。所以随机5点的hough变换存在很严重的系统误差,当取样点分散在椭圆上下左右时,这种误差会比较小,当集中在某个区域时,误差就会非常大。

   

        解决办法就是多采几个点,然后用最小二乘法求解。

下图是10点随机采样的结果:

下图是将所有点一起计算的结果:

可见,采样点越多,拟合度越好。但是一次取样的点越多,这些点落入相同椭圆的概率就越小。这就需要一些手段把椭圆的边缘从噪声中提取出来。至于最小二乘法的椭圆检测算法我将另开一贴来讨论。

 

 

 


 

发布了5 篇原创文章 · 获赞 6 · 访问量 1万+
展开阅读全文

没有更多推荐了,返回首页

©️2019 CSDN 皮肤主题: 大白 设计师: CSDN官方博客

分享到微信朋友圈

×

扫一扫,手机浏览