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

您所在的位置:网站首页 catia椭圆短轴怎么约束 图像处理中的椭圆拟合(一)

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

2023-05-22 12:10| 来源: 网络整理| 查看: 265

原文链接: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