Center detection of circular marking points for opencv camera calibration

1. Background introduction

When scanning the same object from multiple angles , Usually, circular mark points are pasted on objects or other fixed supports to assist in splicing , By calculating the three-dimensional coordinates of the mark points at two angles , Establish corresponding point relationship , utilize SVD Solve the rotation translation matrix . stay Three dimensional reconstruction Before , It is necessary to extract the center coordinates of the mark points on the two-dimensional image , This paper introduces a method of extracting rough contour by using gradient , Then sub-pixel acquisition is performed , Connection profile , The method of ellipse fitting to get the center point . The subpixel contour extraction method described in this paper comes from references [1], About the detailed explanation of the paper, it is also described in other blog posts , This paper focuses on the use of this method to achieve the whole process of mark point center extraction .

2. Center point extraction process

2.1 Rough positioning of target points

The circular sign points with black background and white foreground are shown in the figure 1 Shown , There will be a strong contrast in the image . therefore , Proper threshold binarization can separate the foreground of the marker , The binarized image is arranged in connected domain , According to the connected domain Pixels Number , The aspect ratio of connected domains, etc , Roughly locate the rectangular range of the mark point , Pictured 2 Shown . The connected domain method used is divided into two steps , The first step is to find and number the clusters of locally connected domains , The second step is to traverse the image again , Merge the equivalent connected domains , May refer to [2].

2.2 Subpixel contour extraction

Within each rectangle , The whole pixel 、 Subpixel edge detection , The literature [1] It's an improved canny Edge extraction method .

2.2.1 Calculate the gradient

Use the gray value of the image , At the image level x Calculate the gradient in the direction Gx, In the vertical direction of the image y Calculate the gradient Gy, And solve the amplitude of the gradient 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 Whole pixel and sub-pixel edge points

The literature [1] Based on the analysis of canny Problems in extracting edge points , In order to avoid in canny The error caused by interpolation calculation when solving the gradient of non integer pixels in the process , One of the proposed improvements is to keep only the horizontal or vertical interpolation of edge points . Finally, the edge detection method in the literature will not have the problem of edge oscillation , Pictured 3 Shown .

according to 2.2.1 Gradient magnitude calculated in modG, Judge whether the point is a maximum point in the horizontal or vertical direction , Such as horizontal direction , Judge the current point mod> Left pixel L, And mod<= Right pixel R. Finally, the interpolation is carried out in the horizontal or vertical direction , Depending on which direction the gradient is dominant , That is, if the horizontal interpolation is selected, you need gx > gy, To select vertical interpolation, you need 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 */
}

Find the maximum value of gradient amplitude , And select the direction of interpolation . According to the current point B And before and after A and C spot , A total of three whole pixels modG fitting Quadratic curve , Find the real sub-pixel maximum point on the conic .


Here we mainly derive the source of this formula . Establish the coordinate system as ： With B The point is the origin ,BC The positive direction is x The axis is the pixel coordinate ,y The axis is the gradient amplitude corresponding to the pixel

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 Edge point connection and ambiguity handling

After sub-pixel points are obtained , These points need to be concatenated to form a complete outline . The literature [1] The edge connection method of is also a major feature , Noise points that may appear during edge connection cause "Y" When there is ambiguity such as word splitting , Select the correct edge direction for connection . This paper mainly puts forward the connection criterion of three-point edge contour ：
1） The edge points have similar gradient directions ; That is to say, the two edge points connected , The included angle of their gradient direction is less than 90°. For example, the current point is A, Now judge B Whether the point needs to be connected ,A and B The gradient vectors of are respectively A and B The gradient vectors of are respectively

	//  According to the original description , The judgment is missing from the source code 
	if ((Gx[from] * Gx[to] + Gy[from] * Gy[to]) <= 0.0)
		return 0.0;

2） The line of edge points can separate the bright area from the dark area ; The gradient direction of the edge points is approximately perpendicular to the contour , So the gradient vector rotates clockwise 90°, Should be approximately parallel to the contour , Use judgment conditions ：

among ,AB Is the direction of the outline , And A Gradient vector to 90° The vector dot product of is greater than 0, That is, judge that the included angle between the two is less than 90°. The code section is

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;
}

If B Points exist simultaneously in the backward link A and C, To judge A->B Score of s And A->C Whose score is lower , obviously C leave B It's closer , So the reciprocal of the distance is negative again , The end result will be better than A Point to B Of s smaller , therefore A Point to B The link will be cut off . Empathy , To judge this problem on the forward link “Y” Ambiguity , Judge who gets the higher score .

2.2.4 canny Double threshold processing

Two thresholds are set th_h and th_l, Represents high gradient intensity threshold and low gradient intensity threshold respectively , Above th_h The edge point of is the starting point , Search whether the gradient amplitude of any edge point is less than... Along the forward link direction and the backward link direction respectively th_l, If there is such a point , The link is cut off .

2.2.5 Arrange the final outline

Because the final requirement of a marker is a closed contour , Therefore, it is necessary to select closed contour points from contour points , The meet

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].

// Fitting deviation 
	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++ Realization

On sub-pixel edge point extraction part , The original text has provided the relevant code , In order to realize the center detection of circular mark points , On this basis, the code is modified , In the project .bmp Is the test image ,Center.txt Output center point image coordinates . In recognition, the length limit of the edge contour can be set , For example, greater than 15 And less than 200, And the ellipse fitting error is limited to 2.0.
csdn: Circle mark center extraction c++ engineering
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;
	}
}
 
 reference 
 
[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]  Image connected domain marker 
 [3]  Plane ellipse fitting