Dynamic Granularity Matrix Space Model Based High Robust Multi-Ellipse Center Extraction Method for Camera Calibration

Ellipse center extraction is the important basis of camera calibration with circle arrays calibration board, which directly affects the measurement accuracy of machine vision system. Aiming at the problem of high noise sensitivity of the most ellipse center extraction methods, a high robust multi-ellipse center extraction method is proposed. The idea is to transform the problem of multi-ellipse center extraction into the problems of ellipse sub-pixel edge segmentation and edge based ellipse center clustering. Firstly, the flowchart of the proposed method is introduced. Then, the theory of fuzzy quotient space is introduced to establish the dynamic granularity matrix space model, which describes image segmentation problem as the jumping and transformation of the image by hierarchical structure. Finally, an adaptive multi-ellipse sub-pixel edge segmentation method and an optimal ellipse center extraction method are proposed based on the dynamic granularity matrix space model. Experimental results show that the proposed method has higher accuracy, robustness and lower camera calibration error.


I. INTRODUCTION
Benefiting from the advantages of non-contact and high precision, machine vision [1]- [6] has been widely used in quality detection, three-dimensional reconstruction, visual detection and other fields. The accuracy of camera calibration significantly affects the precision of stereo vision measurement system [7]- [10], therefore, accurate camera calibration needs to be studied for reducing the measurement errors and obtaining the reliable three-dimensional (3D) deformation fields of the object [10]. The task of camera calibration is to determine the relationship between the twodimensional (2D) coordinates of image points and the 3D world coordinates of the object, which takes into account the optical geometry of cameras and the positional relationship between cameras. In order to improve the robustness and calibration accuracy, various calibration boards are used for camera calibration, such as planar chessboard [11], Charuco The associate editor coordinating the review of this manuscript and approving it for publication was Amir Masoud Rahmani . board [12], circular arrays calibration board [13], [14], and so on. In practice, the extraction accuracy of corner points with planar chessboard is usually affected by the lighting environment and noise, leading to obvious fluctuation of calibration results [14]. Although Charuco board can achieve higher calibration accuracy than planar chessboard, Charuco board has higher production complexity and is also affected by extraction accuracy of corner points under uneven illumination. Compared with planar chessboard and Charuco board, due to the advantages of easy detection and high positioning accuracy, the circular arrays calibration board is widely used for robust camera calibration by taking the circular target centers as the matching points of binocular images. However, due to the perspective projection, the circular target will present an ellipse [15] on the binocular image. Therefore, the accuracy of ellipse center extraction directly affects the calibration accuracy, which is very difficult due to the presence of noise, clutter edges, poor lighting conditions, shadow conditions, etc. Ellipse detection is a fundamental technique in image processing field and plays an indispensable role in shape detection and geometric measurement [16]- [20], while the critical issue of detecting multi-ellipse accurately and efficiently is still a challenge work.
Over the past years, many ellipse detection algorithms sprang up and were studied broadly to solve the ellipse detection problem, which can be briefly grouped into several categories: Hough transform method [21]- [24], least-square fitting method [25]- [29], clustering method [30]- [34]. Hough transform has been widely used for detecting geometric primitives such as line segment, circle and ellipse. The basic principle of Hough transform method is to vote arbitrary edge pixels into five-dimensional parameter space, where the detection problem of the given curve in original image is transformed into the problem of finding the local peak in parameter space. Considering the images with strong noise, Lu and Tan [22] proposed an iterative randomized Hough transform method to detect the ellipses. Chen et al. [23] developed a novel ellipse detection method based on Hough transform and edge following method to save the memory and computation time. Combined with artificial neural network, Serra et al. [24] also presented a modified Randomized Hough Transform method to make multiple scans over the acquired images. However, Hough Transform methods suffer from the following legacy problems [22], on one hand, the detected ellipses are often smaller than the real ellipse and sensitive to noise; on the other hand, due to the random sampling, if the selected points are relatively close, it will cause a large deviation of the fitted ellipse. Also, Hough transform requires much effort to tune the model parameters, which will cause heavy computation burden and excessive consumption of memory. Zhang et al. [25] employed least-square ellipse fitting method to estimate the center of the section profile of cylindrical part, which can obtain the cylinder pose precisely. Liang et al. [26] presented a robust ellipse fitting approach by using the L p -norm in the direct least square fitting method. Dong et al. [27] combined the polynomial interpolation algorithm with least-square ellipse fitting method to realize the rapid sub-pixel centroid extraction from a noise-containing spot image. Considering the multiple ellipses detection problem, Grbić et al. [28] proposed a method based on direct least squares method and the random sample consensus algorithm, which performed significantly better for ellipses with noisy edges. The least-square ellipse fitting method has high fitting accuracy, but only one ellipse can be fitted at a time, that is, the image information should be classified and separated before. The general idea of ellipse detection with clustering method is that the pixels are extracted from the arc, then filtered and clustered, and finally the ellipse is fitted by the least-square method [22], [30]- [34]. This method can effectively deal with complex situations such as multiple ellipses, mutual occlusion of ellipses and partial defect of ellipses, which has attracted wide attention. Considering the multiple ellipse fitting problem, Marošević and Scitovski [30] and Scitovski and Marošević [31] proposed a center-based clustering method based on K-means algorithm. Then, Scitovski and Sabo [32] presented a modification of the K-means algorithm for Mahalanobis-circle centers. The centers were determined by iterations of the direct global optimization algorithm, which performed well on the ellipses with clear edges and noisy edges. Also, Li [33] proposed a framework with five processes for the above problem. Singh et al. [34] applied clustering algorithm for image segmentation to obtain the edge pixel of lip, and then the lip ellipse contour was extracted by least-square fitting. It is worth mentioned that in practical applications, ellipse fitting based on clustering method often follows an edge detection step, and thus the effect of edge point errors on the fitting performance cannot be ignored. The K-means method also need set the clustering number and clustering center in practical applications. When there are many ellipses in one image, the application of K-means method will become very difficult.
However, the above ellipse fitting methods [21]- [34] are sensitive to noise, which will reduce the accuracy of camera calibration with circular feature target. In order to achieve higher accuracy of ellipse center extraction and reduce the calibration error, we develop a high robust multi-ellipse center extraction (HRMCE) method based on the dynamic granularity matrix space(DGMS) model. The idea is to transform the problem of multi-ellipse center extraction into the mixed problems of ellipse sub-pixel edge segmentation and edge based ellipse center clustering. For the problem of ellipse sub-pixel edge segmentation, an adaptive sub-pixel edge segmentation method for multi-ellipse is proposed. The region of interest (ROI) of each ellipse in calibration images is extracted by Otsu operator, and the optimal clustering of each ROI can be also obtained based on the DGMS model, which is computed as the adaptive segmentation thresholds for Canny-Zernike operator to obtain the sub-pixel edge points of each ellipse. For the problem of edge based ellipse center clustering, the iterative least-square ellipse fitting method is used to calculate the center points for each ellipse. The DGMS model is used to obtain the optimal clustering of ellipse centers, and extract the optimal ellipse center of each ellipse. The simulation ellipse fitting results under different noises show the effectiveness and higher precision and higher robustness with the proposed HRMCE method. Experimental results of camera calibration demonstrate that compared with other methods, the calibration error with the HRMCE method is significantly reduced. Section 2 explains the basic progress of camera calibration. Section 3 details the principle of the proposed method. Section 4 shows the experimental performance of the proposed method. Sec. 5 summarizes the paper.

II. CAMERA CALIBRATION
Binocular camera calibration is used to characterize the relationship between the actual point of the object and the pixel points of images, which can be presented by camera parameters including camera internal parameters and camera external parameters [1], [2], [9]. Among various famous calibration methods, Zhang's calibration method [9] is widely used in camera calibration. The schematic diagram VOLUME 8, 2020 of binocular vision is shown in Figure 1. O W , X W , Y W , Z W represents the world coordinate system. O l , X l , Y l , Z l and u l , v l denote the left camera coordinate system and its pixel coordinate system, respectively. (O r , X r , Y r , Z r ) and (u r , v r ) denote the right camera coordinate system and its pixel coordinate system, respectively. Point P is any point of the object, o l and o r are the light centers of the left camera and right camera, respectively. P l is the projection point of the point P on the left camera plane. P r is the projection point of the point P on the right camera plane. Then, the line of the light center o l and the point P l must intersect at a point with the line of the light center o r and point P r . P l and P r are oneto-one correspondence. Generally, the calibration target is at the plane z W = 0. According to the ideal pinhole model of camera [9], the transformation relationship of point P from the world coordinate to the pixel coordinate in the image of left camera can be expressed as, where d x , d y represent the size factor of the pixel along the u axis and v axis respectively, (u 0 , v 0 ) represents the main point coordinate of the left camera image, f is the focal length of the camera, f x = f d x , f y = f d y . f x , f y , u 0 , v 0 are only related to the internal structure of the camera. R l T l denotes the camera external parameters, completely determined by the camera coordinate system and the world coordinate system. R l is a rotation matrix with 3 × 3, describing the rotation attitude from the world coordinate to the camera coordinate. T l is a translation vector with 3 × 1, which describes the translation position from the world coordinate to the camera coordinate system.
Let R l = r 1 r 2 r 3 , T l = t l . Equation 1 can be rewritten as, where, A l is the internal parameter matrix of the left camera, H l is the projection matrix of the left camera, which can be written by column vector, Above is a model for ideal lens imaging, but in practical application, lens distortion has a great impact on imaging. Therefore, we introduce a lens imaging model with distortion [9]. If the ideal image point on the normalized image plane is x p , y p , the position of the distortion point is (x d , y d ), and r 2 = x 2 + y 2 ,the position of the ideal image point can be obtained through the following conversion, where (k 1 , k 2 ) is the radial distortion coefficient and (p 1 , p 2 ) is the tangential distortion coefficient.
Since r 1 and r 2 are orthogonal, considering the Equation 3-4, the constraint conditions of internal parameters are obtained, We can see that B is a symmetric matrix. Let vector we can get, By combining the Equations 4-7, we can establish the homogeneous equation Through Equation 8, the camera internal parameters are obtained. Then, the external parameters and the translation vector can be calculated, Similarly, the internal parameters A r and external parameters R r T r of the right camera also can be obtained according to the above steps. According to the above, camera calibration is to obtain the camera internal parameters, f x , f y , u 0 , v 0 , k 1 , k 2 , p 1 , p 2 , the camera external parameters, R l T l . Due to the high robustness, the circular array target is widely used in the camera calibration. However, due to the influence of perspective projection, the circular pattern will present ellipse on camera image, which increases the complexity of feature point extraction. Therefore, the accurate extraction of the ellipse center is very important for the binocular calibration.

III. PRINCIPLE
The idea of the proposed HRMCE method is to transform the center extraction of multi-ellipse into ellipse sub-pixel edge segmentation and ellipse center clustering. Firstly, the dynamic granularity matrix space model is introduced. Then, an adaptive sub-pixel edge segmentation method and optimal ellipse center extraction method are respectively proposed for the above two problems.

A. FLOWCHART OF THE PROPOSED HRMCE METHOD
The flowchart of the proposed HRMCE method is shown in Figure 2. The steps of the HRMCE method are as follows: (1) Obtaining the camera calibration images with circle arrays calibration board.
(2) Using Otsu operator to get the ROI of each ellipse.
(3) By introducing the dynamic granularity matrix space model, the optimal initial population of each ROI can be obtained as the input of the quantum-inspired group leader optimization(QGLO) algorithm [37], [41] to calculate the segmentation thresholds of each ROI. The sub-pixel edge points of each ellipse ROI are calculated by the Canny-Zernike operator.
(4) A number of ellipse center point clusters are obtained by the iterative least-square method with the above sub-pixel edge points.
(5) Based on the DGMS model, the center points of each ellipse are dynamically clustered to get the optimal granularity layer. Compared with the ellipse center with least-square fitting, the one with the smallest distance with the centers is taken as the optimal robust ellipse center. (6) According to the most widely used Zhang's plane calibration method [9], the camera is calibrated with the center point of the circular target as the feature point, and the internal parameters and external parameters of the camera can be obtained.
The proposed HRMCE method has the following advantages. Based on the DGMS model, the problem of multi-ellipse center extraction is transformed into the problems of sub-pixel edge threshold segmentation and ellipse center granularity clustering under different granularities. DGMS model plays a different role in threshold segmentation of ellipse edge segmentation and granularity clustering of ellipse centers. In the edge threshold segmentation, the diverse initial population is obtained based on DGMS model, which provides the initial solution set for the QGLO algorithm to solve the optimal segmentation thresholds.

B. DYNAMIC GRANULARITY MATRIX SPACE MODEL
In this section, we introduce the theory of fuzzy quotient space [35]- [39], and propose the dynamic granularity matrix space model. Considering the consistence of granularity division and image segmentation theory, a dynamic granularity matrix space model is proposed, and the problem of image segmentation can be reconstructed based on the DGMS model. Different from the traditional methods such VOLUME 8, 2020 as threshold segmentation and edge segmentation, the DGMS model describes image segmentation problem as the jumping and transformation of the image at different granularities by hierarchical fuzzy clustering structure [39]. The DGMS model can be used not only for image sub-pixel edge segmentation, but also for ellipse center extraction, which improves the accuracy and robustness of multi-ellipse center extraction.
Definition 1: A triple (X, F, T ) is used to describe the problem. Given an equivalent relation R on X, corresponds to a quotient set [X], thus, the study of the problem (X, F, T ) is transformed into the study of the problem (

[X], [F], [T]). ([X], [F], [T ]) called the quotient space of (X, F, T).
X denotes the domain of the problem, F denotes the attribute function of the domain, T denotes the structural relationship of domain, that is, the relationship among the elements in the domain X . Analyzing or solving problems (X , F, T ) refers to the analysis and study of domain X and its related structures and attributes. Granularity refers to the degree of domain division. Quotient space theory is proposed to analyze and study problems in different granularity spaces, and then synthesize the solutions of the original problems.

Theorem 2: If the problem [A] → [B] has no solution on ([X], [F], [T]), then the problem A → B must have no solution on (X, F, T). Or, if problem A → B has solution on (X, F, T), then the problem [A] → [B] must have solution on quotient space ([X], [F], [T]). Theorem 3: If the problem [A] → [B] has solution on ([X], [F], [T]), let p : X → [X ] be a natural projection,
) is a connected set on X, then the problem A → B must have solution on (X, F, T).
Through Theorem 2, by selecting the proper quotient space and deleting the part without solution, the solution speed can be accelerated. Through Theorem 3, for the proper quotient space, the solution of the original problem can be simplified to the solution of a certain quotient space, and the scale of quotient space is smaller than the original space, which can greatly reduce the computational complexity.
Definition 4: X × X is the product space of X and X , T (X × X ) is the set of fuzzy subsets on the product space, letR ∈ T (X × X ), if satisfies the following conditions, then,R is a fuzzy equivalence relation on X. Theorem 5: LetR be a fuzzy equivalence relation on X , From Definition 1 and Theorem 5, equivalence relation R is only a special case of fuzzy equivalence relatioñ R, that is,R = 0, 1. By introducing the fuzzy equivalence relation, the relation between granules is no longer a definite relation, but a relation under one equivalence degree. Based on the above theory of fuzzy quotient space, we introduce the dynamic granularity matrix space model as follows.  [37]- [38].
In Definition 6, U denotes the studied domain; A denotes the attribute set of granules; V is the set of attribute values of granules, ∃∀a ∈ A, f (a) ∈ V ; I and E respectively reflects the internal attribute and external attribute of the interaction between granules under the same granularity; M represents the functions of the connotation extension operator and the extension connotation operator in the form of granularity matrix.R λ denotes the fuzzy dynamic granularity, reflecting the dynamic characteristics of granules under fuzzy equivalence relation.
Theorem 8: Equivalent relation R λ 1 and equivalent relation R λ 2 respectively correspond to a DGMS model.
Theorem 9: Let {U (λ) |λ ∈ [0, 1] } be a hierarchical structure on U, if given a fuzzy equivalent relationR λ on U, cut-off matrix R λ corresponds to U (λ). The hierarchical structure can be represented by the set of fuzzy equivalent granule Theorem 10: For the dynamic granularity matrix space model U , A, V , (I , E, M ) ,R λ , the following conclusions are equivalent: given a fuzzy equivalent relation; given a hierarchical structure; given a distance function.
According to Theorem 8, through the inclusion relation of quotient sets, the cluster of quotient sets will generate an ordered chain, that is, a hierarchical structure [35]- [38]. According to Theorem 8-9, a fuzzy equivalence relation corresponds to a hierarchical structure, also corresponds to a distance function. For a given distance function, a fuzzy equivalence relation is also defined. Through hierarchical structure, the same problem can obtain the different granularities (or equivalent degrees) and multi-layer clustering structure. Any granularity layer can be used to represent clustering results. When a granular is used to represent the problem space, the granularity is the coarsest; with the granularity refinement, the problem shows more detailed features, which can better describe the problem. The whole information expressed in single layer or single granularity is transformed into attribute information expressed in multi-level and multi-granularity. Therefore, by changing the granularity dynamically, clustering results on different granularity layers can be obtained.

C. ADAPTIVE SUB-PIXEL EDGE SEGMENTATION METHOD FOR MULTI-ELLIPSE
The circle arrays calibration board usually contains dozens of circle targets. The traditional circle target segmentation methods [9] usually use the unified segmentation threshold to extract the edge. However, the unified segmentation threshold will lead to the reduction of segmentation accuracy. For this reason, we propose an adaptive sub-pixel edge segmentation method for multi-ellipse. Firstly, Otsu operator is used to compute the region of interest(ROI) of each circular target. Then, the adaptive segmentation thresholds of each ROI are calculated by DGMS model as the thresholds of Canny-Zernike operator [40]. The ellipse sub-pixel edge of each ROI can be obtained.
Definition 11: Let ROI m n , n = 1, 2, . . . , N n , m = 1, 2, . . . , N m be the mth ellipse ROI in nth calibration image. N n is the number of calibration image, N m is the number of ellipse ROI in each calibration image.
According to the proposed dynamic granularity matrix space model, each ROI can be described by the five tuples X , A, V , (I , E, M ) ,R λ . X is the granular set for each ROI. A = {α, β} denotes the pixel attributes, including gray attribute α and neighborhood gray attribute β. V denotes the set of attribute values, i.e. ∃∀α, β ∈ A, f (α, β) ∈ V . According to pixel attributes A, each ROI can be divided into several fuzzy equivalent clusters, that is, granules. Let connotation of granule, that is, after segmentation, the consistent properties of the region pixels; E is a set of regions with obvious characteristics formed by the consistent pixels. M denotes the relationship between pixels and partitioned regions.
According to Theorem 10, given a distance function, a fuzzy equivalence relation and a hierarchical structure on each ROI can be defined. We define the distance between pixels by normalized distance as follows, where g and g N represent the pixel gray and neighborhood gray respectively. By using ellipse ROI, the global threshold of the ellipse edge is transformed into multiple local thresholds of single ellipse. By introducing DGMS model, the hierarchical structure of each ROI is realized and the optimal clustering also can be obtained. We can see that the dynamic granularity matrix space model transforms the single-layer and single granularity image segmentation problems into multilayer and multi-granularity problems, realizing the dynamic fuzzy image segmentation, that is, the jumping and transformation on different granularity and different granularity layers. By changing the granularity, when the segmentation accuracy is low, the image is described on the coarser granularity layer; when the segmentation accuracy is high, it is necessary to refine the granularity to describe the image on the finer granularity. From hierarchical structure, we can have different clustering results with various granularity layers [40]. Each granularity layer corresponds to one kind of image segmentation results. However, not all of the granularity layers are meaningful for image segmentation. Information entropy [37] is used to measure the effective granularity layer from the hierarchical structure, which can be defined as, where N t denotes the tth information granule in the same granularity layer and |N t | is the number of pixels in the information granule; k denotes the number of information granule; K is the total number of pixels for each ROI. Information gain [37] is used to describe the difference of information entropy between two adjacent granularity layers, shown as, Equation 12 refers to the needed information amount from the coarser granularity layer λ i to the finer granularity layer λ i+1 . The granularity layer has greater information gain means that the granularity layer is more meaningful for image segmentation. Therefore, according to the information entropy method, the optimal granularity layer, that is, optimal clustering results of each ROI can be selected out from the hierarchical structure.
In order to solve the adaptive segmentation thresholds for each ROI, the above optimal clustering results are input into the quantum-inspired group leader optimization(QGLO) algorithm [37], [41] as the initial population, and 2D Otsu is selected as the objective function. Denote the segmentation threshold of pixel gray-scale by s * , and the segmentation threshold of pixel neighborhood gray-scale by t * . Let T upper = ε (s * + t * ) be the upper threshold, ε is a scale factor; T lower = T upper 4 be the lower threshold. We use the adaptive threshold T upper , T lower to improve the Canny-Zernike operator [40] for obtaining the sub-pixel edge points of each ROI, where κ is the number of edge points of the mth ellipse ROI in nth calibration image.

D. OPTIMAL ELLIPSE CENTER EXTRACTION METHOD
In section B, the accurate sub-pixel edge points of each ellipse are obtained. This section introduces the optimal ellipse center extraction method to compute the ellipse center. Figure 3 shows the ellipse diagram. Ellipse can be described by five parameters [25]- [29] including the center point (x 0 , y 0 ), semi-axe (a, b), and major axis deflection angle θ. VOLUME 8, 2020 The general ellipse equation can be presented as, For the original κ ≥ 5 edge points, according to the general ellipse equation and the principle of least-square method, the minimum value of objective function, can be used to determine the parameters (E 1 , Let the partial derivatives of F to each parameter be zero, that is, Thus, the following equations are obtained, x j y j By solving the linear equations, the parameters (E 1 , E 2 , E 3 , E 4 , E 5 ) can be obtained. By substituting the second equation, the parameters of the fitted ellipse can be obtained. The coordinates of the center of the ellipse can be calculated as, From the above, if the least-square ellipse fitting is adopted, at least five different points are needed. Assuming that k points (5 ≤ k ≤ κ) are randomly selected each time. By setting the number of iterations, K , we can have a center point cluster for each ellipse by iterative least-square ellipse fitting method. Let center points of the mth ellipse in nth calibration image be Center m,n = Center k m,n = x0 k m,n , y0 k m,n , k = 1, 2, · · · , K The DGMS model is introduced to fuzzy cluster the ellipse center points. By changing the size of granularity, the cluster structure can be transformed, so that the problem of ellipse center point extraction can be transformed into the transformation and jump of granularity at different levels.
Definition 12: Let Center m,n be the domain, and using a three triples (X c , A c , V c )R λc represents a dynamic granularity matrix space model.
In Definition 12, X c is the granular set formed by all ellipse center points Center m,n . A c denotes the Euclidean distance attribute between center points. V c is the set of attribute values, f A c ∈ V c . According to attribute A c , all ellipse center points of each ROI can be expressed by many fuzzy equivalent clusters represented by granules. The union of all granules can describe the domain Center m,n . When the whole center points are regarded as one granule, the corresponding granularity becomes the finest. While clustering the center points, the granularity becomes coarse. Thus, the problem of optimal ellipse center point extraction can be regarded as the changing process from coarse granularity to fine granularity. According to Theorem 10, given the Euclidean distance of points as d c = abs A c (x i , y i ) − A c x j , y j , we can have a fuzzy equivalence relation and a hierarchical structure about ellipse center points. Given the Euclidean distance, then the fuzzy similarity matrix with symmetry and reflexivity can be expressed as, The effective granularity layer can better reflect the clustering effect of granules, that is, the internal distance between granules is as small as possible, and the external distance between granules is as large as possible. Therefore, the sum of the distances in granules is defined as the density in granules, and the sum of the distances between grains is defined as the separation between granules. The optimal granularity layer has the maximum density within granule and the separation between granules. The criteria for determining the optimal granularity layer can be written as follows, where, Sim avg denotes the density within a granule under granularity λ, Dis avg denotes the separation between granules under granularity λ. Assuming that there are m granules under the granularity, and C λ i represents the ith granule, Let d (i, j) represent the distance between granules, which can be calculated according to Euclidean distance. The average distance between granules is written as, Through Equation 21, we can obtain the optimal clustering, in which the clustering results often contain multiple clusters, that is, there has different cluster centers, S, in the optimal clustering. Different cluster centers represent different ellipse center results, that is, can be represented by EC m,n = EC s m,n = x0 s m,n , y0 s m,n , s = 1, 2, · · · , S . (24) How to select the optimal center point from multiple center points clusters to represent the ellipse center, we present a decision method. Considering that the least-square ellipse fitting method is the most widely used in practical application, and has the characteristics of high precision and high efficiency. Therefore, the ellipse center obtained by the least-square ellipse fitting method is used as a reference value, EC , and all the cluster centers obtained by Equation 21 are compared with the reference value. The cluster center with the least Euclidean distance between them is used as the optimal ellipse center. Let EC o m,n be the optimal ellipse center.
By using the Equation 25, on the one hand, the coarse values can be quickly eliminated to avoid large errors; on the other hand, the extracted optimal ellipse center point can be set near the reference value to improve the accuracy.

IV. EXPERIMENTS AND DISCUSSIONS
In this section, all simulations listed here are implemented in Matlab R2018b on a laptop equipped with 2.50 GHz CPU and 4G RAM memory. Figure 4 shows the binocular vision system with two industrial cameras (model: IDS UI-3370CP-M-GL). Each camera has a resolution of 2048 × 2048 pixels at a frame rate up to 80 fps.
In order to better verify the proposed method, standard ellipse points are generated by computer simulation, and then different noise intensities are added to them to simulate the real situation. As shown in Figure 5, in order to better illustrate the robustness of the proposed method, 4 kinds of noises with different intensities were set in advance. From Figure 5, the ellipse points deviate from the standard ellipse points under different noise intensities. With the increase of noise intensity, the deviation values of ellipse points also become larger. Figure 6 shows the ellipse centers cluster with least-square ellipse fitting method under different numbers of points, using 400 iterations. It can be seen that the more points are selected, the more accurate the center extraction is, while increasing computing burden. Also, the traditional leastsquare fitting with all points is only a special case of the proposed method. The number of selected edge points is set 30 in this experiment.
Then, we investigate the influence of the number of iterations on the calculation of the least-square fitting method. VOLUME 8, 2020   In order to better illustrate the center extraction effect of the simulated ellipse, the HRMCE method is compared with the Hough Transform(HT) method [21]- [24], Least-square ellipse fitting (LS) method [25]- [29], Minimum enclosing circle(MEC) method [42]- [43], K-means method [31]- [34], as shown in Table 1-5, respectively. Error distance can be computed by the Euclidean distance between the real ellipse center and the fitted ellipse center.
From Table 1-5, the HRMCE method achieves the lowest error distance under different noises. Compared with the HT method, the fitting errors with HRMCE method are reduced      HT method, MEC method and K-means method, due to the obvious influence of noise, the errors changes do not increase strictly. For Noise 1, the fitting errors with various methods present that the HRMCE method < HT method < LS method < K-means method < MEC method. For Noise 2, the fitting errors with various methods present that the HRMCE method < K-means method < HT method < LS method < MEC method. Under Noise 3, the fitting errors with various methods present that the HRMCE method < HT method < K-means method < LS method < MEC method. For Noise 4, the fitting errors with various methods present  that the HRMCE method < MEC method < HT method < K-means method < LS method. From Figure 8, the center of the MEC method is determined by the circle composed of the farthest two points. From Table 3, the fitting center error with MEC method is greatly affected by noise. As can be seen from the Figure 8, K-means method divides the points cluster into three categories, while HRMCE method divides the point clusters into four categories. Compared with Kmeans method, HRMCE method does not need to specify the number and center of clustering in advance. From Table 1-5, burden time with various methods present that the HRMCE method > K-means method > MEC method > HT method > LS method. HRMCE method, K-means method and MEC method need more time than LS method to generate the cluster of ellipse center points. HT method requires much more effort to tune the model parameters than LS method. Compared with the MEC method, K-means method requires cluster operation for a large number of ellipse center points. Due to the granularity division and calculation, the HRMCE method has the most burden time, but the ellipse center results obtained are more accurate and more reliable.
The center extraction experiment of the simulated ellipse above shows the high accuracy, effectiveness and robustness of the proposed method. Next, we use the proposed method for binocular calibration experiments and evaluate the influence on the re-projection error. The re-projection error is a geometric error between the projected point and the measured point, which can be represented by the Euclidean distance between them. The re-projection error is chosen as the criterion of binocular calibration. The calibration target with circular array is shown in Figure 9. The distribution of the circle feature in the calibration target is an array with 11 × 9. Five big circles determine the world coordinate system of the target, and the physical distance between the adjacent circle centers is 30 mm. The physical coordinates of 99 centers on the target are used as the feature points of the calibration target. Figure 10-11 show the calibration images of the two cameras. Figure 12 shows the edge segmentation and center extraction results of 3rd calibration image of left camera with the proposed HRMCE method.
In order to evaluate the calibration accuracy, let (x i , y i ) be the world coordinate corresponding to the extracted target center coordinate, (x * i , y * i ) be the known world coordinate of the mark point in the calibration plate, N u be the number of  circular targets in one calibration image, the average error of Euclidean distance, Ave_e, maximum error, Max_e, standard deviation, STD_e, are selected as the evaluation indexes, which can be expressed as,  Table 6 shows the calibration results of left and right cameras with the proposed HRMCE method. In order to better illustrate the calibration accuracy, the HRMCE method is compared with the least-square method [13], True coordinate computation of circle center (TCC) method [44], as shown in Table 7-10, respectively. Table 7-10 are the calibration accuracy analysis with different number of calibration images for left camera.
From Table 7, for 9 calibration images, the mean Ave_e with least square method is 0.0295mm, the mean Max_e is 0.0602mm, and the mean STD_e is 0.0137mm; the mean Ave_e with the proposed HRMCE method is 0.0280mm, the mean Max_e is 0.0565mm, and the mean STD_e is 0.0126mm, reducing by 5.1%, 6.1%, and 8.0%, respectively. Compared with the TCC method, the mean error values with HRMCE method are reduced by 3.1%, 3.4%, and 4.5%, respectively. From Table 8, for 10 calibration images, the mean Ave_e with least square method is 0.0289mm, mean Max_e is 0.0579mm, and mean STD_e is 0.0132mm; the mean error values with the proposed HRMCE method is 0.0279mm, the mean Max_e is 0.0553mm, and the mean STD_e is 0.0120mm, reducing by 3.5%, 4.5%, and 9.1%, respectively. Compared with TCC method, the mean error values with HRMCE method are reduced by 2.7%, 5.3%, and 6.3%, respectively. From Table 9, for 11 calibration images, the mean Ave_e with least square method is 0.0288mm, the mean Max_e is 0.0575mm, and the mean STD_e is 0.0132mm; the mean Ave_e with the proposed HRMCE method is 0.0279mm, the mean Max_e is 0.0550mm, and the mean STD_e is 0.0120mm, reducing by 3.1%, 4.3%, and 9.1% respectively. Compared with the TCC method, the mean error values with HRMCE method are reduced by 3.5%, 5.2%, and 5.5%, respectively. From Table 10, for 12 calibration images, the mean error values with least square method is 0.0288mm, the mean Max_e is 0.0575mm, and the mean STD_e is 0.0130mm; the mean error values with HRMCE method is 0.0276mm, the mean Max_e is 0.0544mm, and the mean STD_e is 0.0118mm, which are reduced by 4.1%, 5.4%, and 9.2% respectively. Compared with the TCC method, the mean error values with HRMCE method are reduced by 2.5%, 4.7%, and 5.6%, respectively. Compared with the least square method and TCC method, the calibration accuracy with the proposed HRMCE method is significantly improved. The least square center extraction method has better fitting effect when there are many acquisition points, but for the calibration target, the available edge pixels are few. Due to the influence of random noise, the center deviation of ellipse will increase. For the TCC method, although the center of ellipse is compensated, it is sensitive to noise, which leads to lower robustness. Especially, for the calibration image with dozens    of circular targets, each circular target is affected by different noises, resulting in larger calibration errors. For the proposed HRMCE method, iterative least square method is used to transform the edge points with less information into a large number of center point clusters, so as to increase the effective information and improve the accuracy and robustness of center extraction. In addition, it can be seen that the least square method is only a special case of the proposed HRMCE method when all edge points participate in the iterative least square. The proposed HRMCE method is a highly robust ellipse center clustering method, which extends the available information and enhances the robustness through iterative method.

V. CONCLUSION
Circular array target is often used for camera calibration. However, due to the perspective projection, the circular target will appear ellipse on calibration image. To improve the accuracy of ellipse center extraction, a high robust multi-ellipse center extraction method for camera calibration is proposed. The idea is to transform the problem of ellipse center extraction into the mixed problem of ellipse sub-pixel edge segmentation and edge based ellipse center clustering. We proposed the dynamic granularity matrix space model based on the theory of fuzzy quotient space. Based on the DGMS model, an adaptive sub-pixel edge segmentation method for multiellipse and an optimal ellipse center extraction method are introduced. Experiments show that the proposed method performs higher precision, higher robustness, and lower calibration error.
In this paper, the center of the ellipse is used as the center of circular target to calibrate the camera, but there is a certain deviation between the two in practice. In future, the deviation will be studied to achieve higher precision camera calibration.