图像处理中的椭圆拟合(opencv)之看的论文的意思的代码

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

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

  1. void FindEllipse(TImage* OrgGrayImg)  
  2. {  
  3. int ImgHeight = OrgGrayImg->nHeight;  
  4.     int ImgWidth = OrgGrayImg->nWidth;  
  5.     unsigned char * Img = OrgGrayImg->pImage; // 输入图像确保是二值图像  
  6.   
  7. srand((unsigned)time(NULL));  
  8.   
  9.     int totalPt = 0;// 用于统计样本点的个数  
  10.   
  11.     for (i = 0; i < ImgHeight - 0; i++)  
  12.     {  
  13.         unsigned char *imgdata = Img + i * ImgWidth;  
  14.         for (j = 0; j < ImgWidth - 0; j++)  
  15.         {  
  16.             if (!imgdata[j])  
  17.                 totalPt ++;  
  18.         }  
  19.     }  
  20.   
  21.     if (totalPt < 5)  
  22.         return;  
  23.   
  24.     POINT * seq;  
  25.     seq = new POINT [totalPt];  
  26.   
  27.     int count = 0;  
  28.     for (i = 0; i < ImgHeight; i++)  
  29.     {  
  30.         unsigned char *data = Img + i * ImgWidth;  
  31.         for (j = 0; j < ImgWidth; j++)  
  32.         {  
  33.             if (!data[j])  
  34.             {  
  35.                 seq[count].x = j;  
  36.                 seq[count].y = i;  
  37.                 count ++;  
  38.             }  
  39.         }  
  40.     }  
  41.   
  42.     double Para[5]; // 存放结果(5个参数A,B,C,D,E)的数组  
  43.     int Angle_V[360]={0};// 椭圆倾斜角参数空间  
  44.     int *Center_XV = new int[ImgWidth];// 椭圆中心点x坐标参数空间  
  45.     int *Center_YV = new int[ImgHeight];// 椭圆中心点y坐标参数空间  
  46.     int *A_axis_V = new int[max(ImgWidth,ImgHeight)/2];// 椭圆长轴参数空间  
  47.     int *B_axis_V = new int[max(ImgWidth,ImgHeight)/2];// 椭圆短轴参数空间  
  48.       
  49.     memset(Center_XV,0,sizeof(int)*ImgWidth);  
  50.     memset(Center_YV,0,sizeof(int)*ImgHeight);  
  51.     memset(A_axis_V,0,sizeof(int)*max(ImgWidth,ImgHeight)/2);  
  52.     memset(B_axis_V,0,sizeof(int)*max(ImgWidth,ImgHeight)/2);  
  53.   
  54.     double Theta,X_c,Y_c,A_axis,B_axis;  
  55.   
  56. int loop = 1;// 成功求出参数的迭代次数  
  57.     int looptop = loop * 1;// 总的迭代次数(也就是控制计算时间的上限,以免陷入无限循环)  
  58. while(loop > 0 && looptop > 0)  
  59. {  
  60. looptop --;  
  61.     int idx;  
  62.     for (count = totalPt; count > 0; count--)// 打乱样本点排列的顺序  
  63.     {  
  64.         POINT ptrtmp;  
  65.         idx = rand() % count;  
  66.           
  67.         ptrtmp = seq[idx];  
  68.         seq[idx] = seq[count-1];  
  69.         seq[count-1] = ptrtmp;        
  70.     }  
  71.   
  72.     double PioMatrix[5*5];  
  73.     for (i = 0; i < 5; i++)  
  74.     {  
  75.         PioMatrix[i*5] = seq[i].x * seq[i].x;  
  76.         PioMatrix[i*5 + 1] = 2 * seq[i].x * seq[i].y;  
  77.         PioMatrix[i*5 + 2] = seq[i].y * seq[i].y;  
  78.         PioMatrix[i*5 + 3] = 2 * seq[i].x;  
  79.         PioMatrix[i*5 + 4] = 2 * seq[i].y;  
  80.     }  
  81.   
  82.     if (GaussJordanInv(PioMatrix,5) == false)// Gauss-Jordan求逆  
  83.         continue;  
  84.     double sum;  
  85.     for (i = 0; i < 5; i++)  
  86.     {  
  87.         sum = 0;  
  88.         for (j = 0; j < 5; j++)  
  89.         {  
  90.             sum +=  -(PioMatrix[i*5 + j]);  
  91.         }  
  92.         Para[i] = sum;  
  93.     }  
  94.   
  95.     if (pow(Para[1],2) - Para[0] * Para[2] > 0)  
  96.             continue;  
  97.   
  98.         if (fabs(Para[0] - Para[2]) < 1e-20)  
  99.             Theta = 1.570796326;  
  100.         else if (Para[0] > Para[2])  
  101.             Theta = 0.5 * (atan(2.0 * Para[1] / (Para[0] - Para[2])) + PI);  
  102.         else  
  103.             Theta = 0.5 * (atan(2.0 * Para[1] / (Para[0] - Para[2])));  
  104.   
  105.         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]);  
  106.         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]);  
  107.         A_axis = 2 * (Para[0] * pow(X_c,2) + Para[2] * pow(Y_c,2) + 2 * Para[1] * X_c * Y_c - 1)   
  108.                         / (Para[0] + Para[2] - sqrt(pow(Para[0] - Para[2],2) + pow(2.0 * Para[1],2)));  
  109.         B_axis = 2 * (Para[0] * pow(X_c,2) + Para[2] * pow(Y_c,2) + 2 * Para[1] * X_c * Y_c - 1)   
  110.                         / (Para[0] + Para[2] + sqrt(pow(Para[0] - Para[2],2) + pow(2.0 * Para[1],2)));  
  111.           
  112.         A_axis = sqrt(A_axis);  //长轴  
  113.         B_axis = sqrt(B_axis);  //短轴  
  114.   
  115.         int AngleTmp = (int)(Theta * 180 / PI + 360 + 0.5) % 360;  
  116.         Angle_V[AngleTmp]++;  
  117.           
  118.         if (X_c < 0 || Y_c < 0 || A_axis < 0 || B_axis < 0)  
  119.             continue;  
  120.         if (X_c >= ImgWidth || Y_c >= ImgHeight || A_axis > max(ImgWidth,ImgHeight)/2 || B_axis > max(ImgWidth,ImgHeight)/2)  
  121.             continue;  
  122.   
  123.         if (X_c >= 0 && X_c < ImgWidth)  
  124.             Center_XV[(int)X_c]++;  
  125.         if (Y_c >= 0 && Y_c < ImgHeight)  
  126.             Center_YV[(int)Y_c]++;  
  127.         if (A_axis >= 0 && A_axis < max(ImgWidth,ImgHeight)/2)  
  128.             A_axis_V[(int)A_axis]++;  
  129.         if (B_axis >= 0 && B_axis < max(ImgWidth,ImgHeight)/2)  
  130.             B_axis_V[(int)B_axis]++;          
  131.         loop--;  
  132. }  
  133.       
  134.     int Angle,Ai,Bi,Cx,Cy;  
  135.     //  Angle  
  136.     int MaxPara = 0;  
  137.     for (i = 0; i < 360; i++)  
  138.     {  
  139.         if (MaxPara < Angle_V[i])  
  140.         {  
  141.             MaxPara = Angle_V[i];  
  142.             Angle = i;  
  143.         }  
  144.     }  
  145.     //  Cy  
  146.     MaxPara = 0;  
  147.     for (i = 0; i < ImgHeight; i++)  
  148.     {  
  149.         if (MaxPara < Center_YV[i])  
  150.         {  
  151.             MaxPara = Center_YV[i];  
  152.             Cy = i;  
  153.         }  
  154.     }  
  155.     //  Cx  
  156.     MaxPara = 0;  
  157.     for (i = 0; i < ImgWidth; i++)  
  158.     {  
  159.         if (MaxPara < Center_XV[i])  
  160.         {  
  161.             MaxPara = Center_XV[i];  
  162.             Cx = i;  
  163.         }  
  164.     }  
  165.     //  Ai  
  166.     MaxPara = 0;  
  167.     for (i = 0; i < max(ImgWidth,ImgHeight)/2; i++)  
  168.     {  
  169.         if (MaxPara < A_axis_V[i])  
  170.         {  
  171.             MaxPara = A_axis_V[i];  
  172.             Ai = i;  
  173.         }  
  174.     }  
  175.     //  Bi  
  176.     MaxPara = 0;  
  177.     for (i = 0; i < max(ImgWidth,ImgHeight)/2; i++)  
  178.     {  
  179.         if (MaxPara < B_axis_V[i])  
  180.         {  
  181.             MaxPara = B_axis_V[i];  
  182.             Bi = i;  
  183.         }  
  184.     }  
  185.   
  186.     delete[] Center_XV;  
  187.     delete[] Center_YV;  
  188.     delete[] A_axis_V;  
  189.     delete[] B_axis_V;  
  190.   
  191.   
  192.     double sma = SinMem[Angle];  
  193.     double cma = CosMem[Angle];  
  194.     for (int n = 0; n < 360; n++)  
  195.     {  
  196.         i = (Bi) * CosMem[360 - n];  
  197.         j = (Ai) * SinMem[360 - n];  
  198.           
  199.         int x,y;  
  200.         x = (j * cma - i * sma) + Cx;  
  201.         y = (i * cma + j * sma) + Cy;  
  202.           
  203.         Mask[y * ImgWidth + x] = 0;  
  204.     }  
  205.     delete[] Mask;  
  206.     delete[] seq;  
  207. }  

测试结果:

原图:

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

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

原图:

拟合结果:

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

   

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

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

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

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

转:http://blog.csdn.net/easecode/article/details/21188657

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