原文链接: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点随机采样的结果：
下图是将所有点一起计算的结果：
可见，采样点越多，拟合度越好。但是一次取样的点越多，这些点落入相同椭圆的概率就越小。这就需要一些手段把椭圆的边缘从噪声中提取出来。至于最小二乘法的椭圆检测算法我将另开一贴来讨论。
