图像处理中的椭圆拟合(一) |
您所在的位置:网站首页 › catia椭圆短轴怎么约束 › 图像处理中的椭圆拟合(一) |
原文链接:https://blog.csdn.net/easecode/article/details/21188657 图像处理中的椭圆检测用处还是挺多的,找到这里来的同学大多是想用椭圆检测来解决某些实际问题吧,所以我就不做介绍,直奔主题。我研究这块也有一段时间了,也查找了挺多资料,貌似通用的椭圆算法还没有,不排除我孤陋寡闻了。前辈提出的算法适用范围比较有限,这个”比较有限“是相对直线检测来说的。但直接用Hough变换来找椭圆几乎是不可能的事,在5维空间里做”投票“,想想都觉得可怕。于是有人想到用随机Hough变换。这是一种很合理的方法,我就是这么做的,不过这种方法有个不足之处,后面会讲到。这里先介绍这方法的流程。 二次曲线的一般方程为:,其中(x,y)为图像的坐标空间,B、C、D、E、F是二次曲线的参数。当满足时二次曲线为椭圆。方程中需要求解的参数有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点随机采样的结果: 下图是将所有点一起计算的结果: 可见,采样点越多,拟合度越好。但是一次取样的点越多,这些点落入相同椭圆的概率就越小。这就需要一些手段把椭圆的边缘从噪声中提取出来。至于最小二乘法的椭圆检测算法我将另开一贴来讨论。 |
今日新闻 |
推荐新闻 |
CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |