Opencv相机标定之圆形标识点中心检测

1. 背景介绍

多个角度扫描同一个物体时，通常会在物体上或者其他固定支架上粘贴圆形标志点来辅助进行拼接，通过计算两个角度下的标志点三维坐标，建立对应点关系，利用SVD求解旋转平移矩阵。在三维重建之前，需要在二维图像上提取标志点中心坐标，本文讲解一种利用梯度提取粗轮廓，然后进行亚像素获取，连接轮廓，椭圆拟合得到中心点的方法。本文所述亚像素轮廓提取方法来自于参考文献[1]，关于论文的详解其他博文中也有讲述，本文着重是利用该方法实现全流程的标志点中心提取。

2. 中心点提取流程

2.1 目标点粗定位

使用黑色背景白色前景的圆形标志点如图1所示，在图像上会表现出强烈的对比反差。因此，恰当的阈值二值化就能将标志点的前景分离出来，将二值化后的图像进行连通域整理，根据连通域占据像素数量，连通域的横纵比等，粗略的定位出标志点所在的矩形范围，如图2所示。所用的连通域方法分为两步，第一步是各个局部连通域团的查找和编号，第二步是再次遍历图像，将等价的连通域合并，可参考[2]。

2.2 亚像素轮廓提取

在每一个矩形范围内，进行整像素、亚像素的边缘提取，文献[1]是一种改进的canny边缘提取方法。

2.2.1 计算梯度

利用图像的灰度值，在图像水平x方向计算梯度Gx，在图像垂直方向y计算梯度Gy，并求解梯度的幅值modG。

void compute_gradient(double * Gx, double * Gy, double * modG, uchar * image, int X, int Y)
{
    
	int x, y;
	if (Gx == NULL || Gy == NULL || modG == NULL || image == NULL)
		error("compute_gradient: invalid input");
	// approximate image gradient using centered differences
	for (x = 1; x<(X - 1); x++)
	for (y = 1; y<(Y - 1); y++)
	{
    
		Gx[x + y*X] = (double)image[(x + 1) + y*X] - (double)image[(x - 1) + y*X];
		Gy[x + y*X] = (double)image[x + (y + 1)*X] - (double)image[x + (y - 1)*X];
		modG[x + y*X] = sqrt(Gx[x + y*X] * Gx[x + y*X] + Gy[x + y*X] * Gy[x + y*X]);
	}
}

2.2.2 整像素和亚像素边缘点

文献[1]在分析canny提取边缘点时存在的问题，为了规避在canny流程中对非整像素点求解梯度时插值计算导致的误差，提出的改进点之一就是保持只在水平或者垂直进行边缘点插值计算。最终文献的边缘提取方法不会出现边缘振荡的问题，如图3所示。

根据2.2.1中计算的梯度幅值modG， 在水平或者垂直方向上判断该点是否是极大值点，比如水平方向，判断当前点mod>左像素点L，且mod<=右像素点R。最终到底在水平还是垂直方向上进行插值，取决于哪个方向的梯度占主导，即选择水平方向插值的话需要gx > gy，选择垂直方向插值的话需要gy > gx。

double mod = modG[x + y*X];     /* modG at pixel */
double L = modG[x - 1 + y*X];   /* modG at pixel on the left */
double R = modG[x + 1 + y*X];   /* modG at pixel on the right */
double U = modG[x + (y + 1)*X]; /* modG at pixel up */
double D = modG[x + (y - 1)*X]; /* modG at pixel below */
double gx = fabs(Gx[x + y*X]);  /* absolute value of Gx */
double gy = fabs(Gy[x + y*X]);  /* absolute value of Gy */
if (greater(mod, L) && !greater(R, mod) && gx >= gy)
{
    
	Dx = 1; /* H */
}
else if (greater(mod, D) && !greater(U, mod) && gx <= gy)
{
    
	Dy = 1; /* V */
}

在找到梯度幅值的极大值点，并且选择插值的方向后。根据当前点B和前后的A和C点，共三个整像素点的modG拟合二次曲线，在二次曲线上去寻找真正的亚像素极大值点。


这里主要推导一下该公式的来源。将坐标系建立为：以B点为原点，BC正向为x轴为像素坐标，y轴为像素对应的梯度幅值

if (Dx > 0 || Dy > 0)
{
    
	/* offset value is in [-0.5, 0.5] */
	double a = modG[x - Dx + (y - Dy) * X];
	double b = modG[x + y     * X];
	double c = modG[x + Dx + (y + Dy) * X];
	double offset = 0.5 * (a - c) / (a - b - b + c);
	/* store edge point */
	Ex[x + y*X] = x + offset * Dx;
	Ey[x + y*X] = y + offset * Dy;
}

2.2.3 边缘点连接以及歧义情况的处理

将亚像素的点获取之后，需要将这些点串联在一起形成一个完整的轮廓。文献[1]的边缘连接方法也是一大特色，针对边缘连接时可能出现的噪声点导致"Y"字分路等歧义情况时，选择正确的边缘方向进行连接。主要提出了三点边缘轮廓的连接准则：
1）边缘点具有相似的梯度方向；也就是说连接起来的两个边缘点，他们的梯度方向的夹角小于90°。比如当前点为A，现在判断B点是否需要连接，A和B的梯度向量分别为A和B的梯度向量分别为

	// 根据原文描述，源代码遗漏了该判断
	if ((Gx[from] * Gx[to] + Gy[from] * Gy[to]) <= 0.0)
		return 0.0;

2）边缘点连线能够将亮区域和暗区域分离开；边缘点的梯度方向近似的垂直于轮廓，那么梯度向量顺时针旋转90°，应该近似与轮廓平行，利用判断条件：

其中，AB是轮廓的方向，与A梯度向量转90°的向量点乘大于0，即判断两者夹角小于90°。代码部分为

dx = Ex[to] - Ex[from];
dy = Ey[to] - Ey[from];
if ((Gy[from] * dx - Gx[from] * dy) * (Gy[to] * dx - Gx[to] * dy) <= 0.0)
	return 0.0; /* incompatible gradient angles, not a valid chaining */

if ((Gy[from] * dx - Gx[from] * dy) >= 0.0)
	return  1.0 / dist(Ex[from], Ey[from], Ex[to], Ey[to]); /* forward chaining */
else
	return -1.0 / dist(Ex[from], Ey[from], Ex[to], Ey[to]); /* backward chaining */

double s = chain(from, to, Ex, Ey, Gx, Gy, X, Y);  /* score from-to */

if (s > fwd_s)	/* a better forward chaining found */
{
    
	fwd_s = s;  /* set the new best forward chaining */
	fwd = to;
}
if (s < bck_s)	/* a better backward chaining found */
{
    
	bck_s = s;  /* set the new best backward chaining */
	bck = to;
}

假如B点的后向链路中同时存在A和C，就要判断A->B的得分s与A->C的得分谁更低，显然C离B的距离更近，因此距离的倒数再取负，最终结果会比A指向B的s更小，因此A指向B的链路会被切断。同理，若要判断正向链路上的这种“Y”歧义情况，判断谁的得分更高。

2.2.4 canny双阈值处理

设置了两个阈值th_h和th_l，分别表示高梯度强度阈值与低梯度强度阈值，以高于th_h的边缘点为起始点，分别沿着前向链路方向和后向链路方向搜索是否有边缘点的梯度幅值小于th_l，若存在这样点，则链路被切断。

2.2.5 整理最终的轮廓

因为标志点最终需要的是首尾封闭的轮廓，因此需要从轮廓点中挑选出封闭的轮廓点，即满足

curve k is closed if x[curve_limits[k]] == x[curve_limits[k+1] - 1] and
y[curve_limits[k]] == y[curve_limits[k+1] - 1].

//拟合偏差
	for (int i = 0; i < nPt; i++)
	{
    
		Px = vEdgePtList[i].dX;
		Py = vEdgePtList[i].dY;
		dFitDev += abs(Px*Px + dEllipsePara[0] * Px*Py + dEllipsePara[1] * Py*Py + dEllipsePara[2] * Px + dEllipsePara[3] * Py + dEllipsePara[4]);
	}
	dFitDev /= nPt;

3. C++实现

关于亚像素边缘点提取的部分，原文已提供了相关的代码，本文为了实现圆形标志点的中心检测，在此基础上修改了代码，工程中的.bmp是测试图像，Center.txt输出的中心点图像坐标。在识别中可以通过设置边缘轮廓的长度限制，比如大于15且小于200，并且限制椭圆拟合误差小于2.0。
csdn: 圆形标志点中心提取c++工程
github: Reference Point Detection

if (vEdgePt.size() > 15 && vEdgePt.size() < 200)
{
    
	double dCx = 0, dCy = 0;
	double dFitDev = 0;
	FitEllipsePara(dCx, dCy, dFitDev, vEdgePt);
	if (dFitDev < 2.0)
	{
    
		dCx += iX0;
		dCy += iY0;
		fout << dCx << " " << dCy << endl;
	}
}
 
参考文献
 
[1] Gioi R G V , Randall G . A Sub-Pixel Edge Detector: an Implementation of the Canny/Devernay Algorithm[J]. Image Processing On Line, 2017, 7:347-372.
 [2] 图像连通域标记
 [3] 平面椭圆拟合