3.1 Edge coarse localisation based on the Canny operator

1. (1) Gaussian noise reduction of grey image

The original grey images are smoothed and denoised by two-dimensional Gaussian function.

(1) \begin{align}G(u,v) = \frac{1}{{2\pi {\sigma ^2}}}\exp \left( - \frac{{{u^2} + {\nu ^2}}}{{2{\sigma ^2}}}\right) \end{align}

where $\sigma $ is the standard deviation of the Gaussian function. The output image after Gaussian smoothing is obtained as follows.

(2) \begin{align}I(u,v) = G(u,v)g(u,v) \end{align}

where $g(u,v)$ is original grey image.

1. (2) Calculation of image gradient amplitude and gradient direction

   The gradient amplitude image $H(u,v)$ and gradient direction $t(u,v)$ are defined as:

   (3) \begin{align}H(u,v) = \sqrt {H_u^2(u,v) + H_v^2(u,v)} \end{align}
   (4) \begin{align}t(u,v) = \arctan \left( {\frac{{{H_v}(u,v)}}{{{H_u}(u,v)}}} \right) \end{align}

2. (3) Image gradient assignment

The pixel point to be processed is selected as the centre, and the interpolation of the two adjacent gradient amplitudes in the current pixel’s gradient direction is calculated. If the magnitude of the gradient is more than the magnitudes of the two interpolations, the point is classified as an edge point and assigned a value of 1. Otherwise, it is assigned a value of 0.

1. (4) Edge determination

The edge is detected and connected based on the high and low double thresholds. The high threshold is ${T_{\max }}$ and the low threshold is ${T_{\min }}$ . If the gradient intensity value of a pixel in the image is greater than ${T_{\max }}$ , the point is a strong edge point. If the gradient intensity value of a pixel is less than ${T_{\min }}$ , the pixel is a non-edge point. If the gradient intensity value of a pixel is between ${T_{\min }}$ and ${T_{\max }}$ , the point is a weak edge point. The weak edge point may be an edge point or a non-edge point caused by noise, which needs to be further determined by combining strong edge points. Following the above procedures, an image of the unit circle with coarse edge localisation is obtained.

3.2 Sub-pixel edge extraction based on improved Zernike moments

1. (1) Zernike moment

Zernike polynomial is defined as:

(5) \begin{align}{V_{nm}}(\rho ,\omega ) = {S_{nm}}\left( \rho \right){e^{im\omega }} \end{align}

where $n$ , $m$ are integers with $n \geqslant 0,{\rm{ }}n \geqslant \left| m \right|$ is even number and $n \geqslant \left| m \right| \gt 0$ . The expression of ${S_{nm}}\left( \rho \right)$ is as follows.

(6) \begin{align} {S_{nm}}\left( {\rho} \right) = \mathop {\mathop \sum \nolimits^{} }\limits_{k = 0}^{{{\left( {n - |m|} \right)} \mathord{\left/ {\vphantom {{\left( {n - |m|} \right)} 2}} \right. } 2}} {{{{\left( { - 1} \right)}^k}\left( {n - k} \right)!{\rho ^{n - 2k}}} \over {k!\left( {{{n + |m|} \over 2} - k} \right)!\left( {{{n - |m|} \over 2} - k} \right)!}} \end{align}

where $\rho $ denotes the distance from the coordinate point $\left( {u,v} \right)$ in the Cartesian coordinate system to the origin, and $\omega $ is the angle between $\rho $ and the $u$ -direction.

The Zernike moments of the image $f\left( {u,v} \right)$ are defined as:

(7) \begin{align}Ze{r_{nm}} = {{n + 1} \over \pi }\int\!\!\!\int\limits_{{u^2} + {v^2} \le 1} {f(u,v)V_{nm}^*(\rho ,\omega )dud} v \end{align}

Under the discrete condition, the Zernike moments can be expressed in the digital image by Equation (8).

(8) \begin{align}Ze{r^*}_{nm} = \sum\limits_u {\sum\limits_v {{f^*}(u,v)V_{nm}^*(\rho ,\omega )} } \end{align}

where $V_{nm}^*$ is the conjugate matrix of Zernike polynomials in complex domain.

If the image is rotated by an angle $\omega $ , the relationship between the Zernike polynomials satisfies $Ze{r'_{nm}} = Ze{r_{nm}}{e^{ - im\omega }}$ . Therefore, when the image rotates, only the phase angle is changed and the phase modulus of the Zernike matrix remains unchanged. This characteristic is the rotational invariance of the Zernike moment. By using a rotated image, edge parameters can be more easily solved, thereby achieving sub-pixel edge localisation.

1. (2) Improved the calculation of Zernike moment

The traditional sub-pixel edge detection based on Zernike moment uses a step model at the edge, but the edge intensity of the actual image doesn’t change stepwise but gradually. Therefore, the cumulative distribution function of the edge model is built.

(9) \begin{align}I\left( {u'} \right) = h + \frac{k}{2}\left[1 + erf \left(\frac{{u' - l}}{{\sqrt 2 \sigma }}\right)\right] \end{align}

In Equation (9), $u'$ is the new coordinate axis after image rotation, which is perpendicular to the edge. $\sigma $ is the standard deviation, which can change the edge width. l is the distance from edge to $u'$ -axis. In order to simplify the calculation process, the slope edge model is used instead. The functional expression of the slope model is as follows.

(10) \begin{align} I\left( {u'} \right) = \left\{ {\begin{array}{ll} h {}& {u' \le d - \varsigma }\\[8pt] {h + g\frac{{u' - (d - \varsigma )}}{{2\varsigma }}} & {}{d - \varsigma \lt u' \le d + \varsigma }\\[8pt] {h + g} & {}{d + \varsigma \lt u'}\end{array}} \right. \end{align}

where $2\varsigma $ is the horizontal width of the slope, $d$ is distance from the edge to the v-axis, and h and h + g represent the grey value of image background and image, respectively. At this point, the edge model changes from a step function model to a ramp function model, as shown in Fig. 3.

Figure 3.

Step-edge model and ramp-edge model.

To facilitate integration, the defining equation of Zernike moments can be written in the form of right-angle coordinates, and the expression of $Ze{r_{nm}}$ in the rotated coordinate system is

(11) \begin{align}Ze{r_{nm}} = \int\!\!\!\int_{{{u'}^2} + {{v'}^2} \le 1} f (u,v)V_{nm}^*(u',v')du'dv' \end{align}

According to the slope model, the calculation of $Ze{r{^\prime}_{nm}} $ for the new model can be divided into three regions, as shown in Fig. 3, and the expressions are

(12) \begin{align}Ze{r'_{nm}} = Ze{r'_{nm1}} + Ze{r'_{nm2}} + Ze{r'_{nm3}} \end{align}

where $Ze{r'_{nm1}}$ , $Ze{r'_{nm2}}$ and $Ze{r'_{nm3}}$ are

(13) \begin{align}\left\{\begin{array}{l} \displaystyle Zer^{\prime}_{nm1} = h \int^{1}_{-1} \int^{\sqrt{1-u^{2}}}_{\sqrt{1-u^{2}}}\bar{V}_{nm}(u^{\prime},v^{\prime})du^{\prime}dv^{\prime}\\[9pt] \displaystyle Zer^{\prime}_{nm2} = \frac{g}{2\varsigma} \int^{d+\varsigma}_{d-\varsigma}(u^{\prime}- d + \varsigma) \int^{\sqrt{1-u^{2}}}_{\sqrt{1-u^{2}}}\bar{V}_{nm}(u^{\prime},v^{\prime})du^{\prime}dv^{\prime}\\[9pt] \displaystyle Zer^{\prime}_{nm3} = g\int_{d+\varsigma}^{1} \int^{\sqrt{1-u^{\prime 2}}}_{-\sqrt{1-u^{\prime 2}}}\bar{V}_{nm}(u^{\prime},v^{\prime})du^{\prime}dv^{\prime} \end{array}\right.\end{align}

When $n = 1$ and $m = 1$ , according to Equation (12), $Ze{r'_{11}}$ is obtained.

(14) \begin{align}\begin{array}{c}Ze{{r'}_{11}} = \frac{g}{{24\varsigma }}\left\{ 3\arcsin\!(d + \varsigma ) - 3\arcsin\!(d - \varsigma ) - [5(d - \varsigma ) - 2{(d - \varsigma )^3}]\sqrt {1 - {{(d - \varsigma )}^2}}\right. \\[5pt] \left.{\rm{ }} + [5(d + \varsigma ) - 2{(d + \varsigma )^3}]\sqrt {1 - {{(d + \varsigma )}^2}} \right\} \end{array} \end{align}

When $n = 2$ and $m = 0$ , $Ze{r'_{20}}$ is obtained.

(15) \begin{align}Ze{r'_{20}} = \frac{g}{{15\varsigma }}\{ {[1 - {(d - \varsigma )^2}]^{5/2}} - {[1 - {(d + \varsigma )^2}]^{5/2}}\} \end{align}

When ω tends to 0, Equations (16) and (17) are obtained as follows.

(16) \begin{align}\mathop {\lim }\limits_{u \to \infty } Ze{r'_{11}} = \frac{{2g}}{3}{(1 - {d^2})^{\frac{3}{2}}} \nonumber \\[-30pt] \end{align}

(17) \begin{align}\mathop {\lim }\limits_{u \to \infty } Ze{r'_{20}} = \frac{{2gd}}{3}{(1 - {d^2})^{\frac{3}{2}}} \end{align}

From Equations (16) and (17), it can be seen that when $\varsigma $ tends to 0, the results are infinitely close to the results in the step model. It can be said that the step model is a special case of the slope model. In the traditional step model algorithm, d is related to the ratio of two Zernike moments, so the relationship between the two in the ramp model can be approximated as:

(18) \begin{align}\frac{{Ze{r_{20}}}}{{Ze{r_{11}}}} \approx d - {\varsigma ^2}l - \frac{1}{4}{\varsigma ^2}{l^2} \end{align}

According to Equation (18), d is derived from Equation (18) as:

(19) \begin{align}d = \frac{2}{{{\varsigma ^2}}}\left[1 - {\varsigma ^2} - \sqrt {{{({\varsigma ^2} - 1)}^2} - \frac{{{\varsigma ^2}Ze{r_{20}}}}{{Ze{{r'}_{11}}}}}\right] \end{align}

When $\varsigma $ tends to 0, the limit of d is shown in Equation (20), and it can be seen from the expression that d is infinitely close to the result in the step model at this time.

(20) \begin{align}\mathop {\lim }\limits_{\varsigma \to \infty } d = \frac{{Ze{r_{20}}}}{{Ze{{r'}_{11}}}} \end{align}

Let Re denote the real part and Im denote the imaginary part. It can be deduced that

(21) \begin{align}\alpha = \arctan \left( {\frac{{{\mathop{\rm Im}\nolimits} \left[ {Ze{r_{11}}} \right]}}{{{\mathop{\rm Re}\nolimits} \left[ {Ze{r_{11}}} \right]}}} \right) \end{align}

According to Equation (21), the edge normal direction of the image $f(u,v)$ can be calculated as:

(22) \begin{align}f(u,v) = \left[ \begin{array}{l}{n_u}\\{n_v}\end{array} \right] = \left[ \begin{array}{l}\cos \alpha \\\sin \alpha \end{array} \right] \end{align}

Using the pixel-level edge localisation point $({u_0},{v_0})$ as the coordinate origin, the final sub-pixel edge coordinates in the actual image can be obtained.

(23) \begin{align}\left[ {\begin{array}{*{20}{c}}{{u_{\rm{e}}}}\\{{v_{\rm{e}}}}\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}{{u_0}}\\{{v_0}}\end{array}} \right] + \frac{{Nd}}{2}\left[ {\begin{array}{*{20}{c}}{\cos \alpha }\\{\sin \alpha }\end{array}} \right] \end{align}

Where N is template size, and N = 7

During the camera calibration, the centres of the solid circles on the calibration board are arranged in a random sequence. To establish the one-to-one correspondence between the three-dimensional coordinates of circle centres and their pixel coordinates in the images, it is necessary to rearrange the pixel coordinates of these solid circle centres. Therefore, a circle centre sorting method based on triangular marker is introduced. As depicted in Fig. 4, there is an isosceles right-angled triangular feature at the upper-left corner of the calibration board. The hypotenuse of this isosceles right-angled triangle serves as a reference from which the closest circle is selected as the first circle, determining the position of circle centre 1. Following this, the positions of the remaining circle centres are sequentially arranged, completing the sorting process for all circle centres.

Figure 4.

Circle centre location and sorting.

4.1 Intrinsic and extrinsic parameters model

According to the multi-cameras detection system built in Fig. 1, the model of multi-cameras coordinate system is built as shown in Fig. 5.

Figure 5.

The coordinate system of multi-cameras.

In Fig. 5, ${o_{c1}}$ , ${o_{c2}}$ , ${o_{c3}}$ respectively represent optical centres of camera 1, camera 2 and camera 3. ${o_{c1}} - {x_{c1}}{y_{c1}}{z_{c1}}$ , ${o_{c2}} - {x_{c2}}{y_{c2}}{z_{c2}}$ , ${o_{c3}} - {x_{c3}}{y_{c3}}{z_{c3}}$ respectively represent the camera coordinate systems of the camera 1, camera 2 and camera 3. ${O_{c1}} - {X_{c1}}{Y_{c1}}$ , ${O_{c2}} - {X_{c2}}{Y_{c2}}$ , ${O_{c3}} - {X_{c3}}{Y_{c3}}$ . respectively represent the image coordinate systems of the camera 1, camera 2 and camera 3. ${o_1} - {u_1}{v_1}$ , ${o_2} - {u_2}{v_2}$ , ${o_3} - {u_3}{v_3}$ represent the pixel coordinate systems of camera 1, camera 2 and camera 3. p is any point in the world coordinate system. The homogeneous coordinates of p are denoted by $p\left( {{x_w},{y_w},{z_w},1} \right)$ . ${P_{c1}}$ , ${P_{c2}}$ , ${P_{c3}}$ are the imaging points of p on the image coordinate system of camera 1, camera 2 and camera 3, respectively.

The homogeneous coordinates of imaging points of p in the camera 1, camera 2 and camera 3 image coordinate systems are denoted by ${p_{c1}}\left( {{x_{c1}},{y_{c1}},{z_{c1}},1} \right)$ , ${p_{c2}}\left( {{x_{c2}},{y_{c2}},{z_{c2}},1} \right)$ , and ${p_{c3}}\left( {{x_{c3}},{y_{c3}},{z_{c3}},1} \right)$ , respectively.

(24) \begin{align}\left\{ \begin{array}{c}{{\boldsymbol{p}}_{c1}}^{\rm{T}} = \left[ {\begin{array}{*{20}{c}}{{x_{c1}}}\\{{y_{c1}}}\\{{z_{c1}}}\\1\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}{_w^{c1}{r_{11}}}\quad {}{_w^{c1}{r_{12}}} {}\quad{{}_w^{c1}{r_{13}}}\quad {}{{}_w^{c1}{\tau _x}}\\[3pt] {{}_w^{c1}{r_{21}}}\quad {}{{}_w^{c1}{r_{22}}}\quad {}{{}_w^{c1}{r_{23}}}\quad {}{{}_w^{c1}{\tau _y}}\\[3pt] {{}_w^{c1}{r_{31}}}\quad {}{{}_w^{c1}{r_{32}}}\quad {}{{}_w^{c1}{r_{33}}}\quad {}{{}_w^{c1}{\tau _z}}\end{array}} \right] \left[ {\begin{array}{*{20}{c}}{{x_w}}\\{{y_w}}\\{{z_w}}\\1\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}{{}_w^{c1}{\boldsymbol{R}}}\quad {}{{}_w^{c1}{\unicode{x1D71E}}}\end{array}} \right] {{\boldsymbol{p}}^{\rm{T}}} = {}_w^{c1}{\boldsymbol{M}}{{\boldsymbol{p}}^{\rm{T}}}\\[25pt] {{\boldsymbol{p}}_{c2}}^{\rm{T}} = \left[ {\begin{array}{*{20}{c}}{{x_{c2}}}\\{{y_{c2}}}\\{{z_{c2}}}\\1\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}{{}_w^{c2}{r_{11}}}\quad {}{{}_w^{c2}{r_{12}}}\quad {}{{}_w^{c2}{r_{13}}}\quad {}{{}_w^{c2}{\tau _x}}\\[3pt] {{}_w^{c2}{r_{21}}}\quad {}{{}_w^{c2}{r_{22}}}\quad {}{{}_w^{c2}{r_{23}}}\quad {}{{}_w^{c2}{\tau _y}}\\[3pt] {{}_w^{c2}{r_{31}}}\quad {}{{}_w^{c2}{r_{32}}}\quad {}{{}_w^{c2}{r_{33}}}\quad {}{{}_w^{c2}{\tau _z}}\end{array}} \right] \left[ {\begin{array}{*{20}{c}}{{x_w}}\\{{y_w}}\\{{z_w}}\\1\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}{{}_w^{c2}{\boldsymbol{R}}}\quad {}{{}_w^{c2}{\unicode{x1D71E}}}\end{array}} \right] {{\boldsymbol{p}}^{\rm{T}}} = {}_w^{c2}{\boldsymbol{M}}{{\boldsymbol{p}}^{\rm{T}}}\\[25pt] {{\boldsymbol{p}}_{c3}}^{\rm{T}} = \left[ {\begin{array}{*{20}{c}}{{x_{c3}}}\\{{y_{c3}}}\\ {{z_{c3}}}\\1\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}{{}_w^{c3}{r_{11}}}\quad {}{{}_w^{c3}{r_{12}}}\quad {}{{}_w^{c3}{r_{13}}}\quad {}{{}_w^{c3}{\tau _x}}\\[3pt] {{}_w^{c3}{r_{21}}}\quad {}{{}_w^{c3}{r_{22}}}\quad {}{{}_w^{c3}{r_{23}}}\quad {}{{}_w^{c3}{\tau _y}}\\[3pt] {{}_w^{c3}{r_{31}}}\quad {}{{}_w^{c3}{r_{32}}}\quad {}{{}_w^{c3}{r_{33}}}\quad {}{{}_w^{c3}{\tau _z}}\end{array}} \right] \left[ {\begin{array}{*{20}{c}}{{x_w}}\\{{y_w}}\\{{z_w}}\\1\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}{{}_w^{c3}{\boldsymbol{R}}} {}{{}_w^{c3}{\unicode{x1D71E}}}\end{array}} \right] {{\boldsymbol{p}}^{\rm{T}}} = {}_w^{c3}{\boldsymbol{M}}{{\boldsymbol{p}}^{\rm{T}}}\end{array} \right. \end{align}

where ${}_w^{c1}{\boldsymbol{R}}$ , ${}_w^{c2}{\boldsymbol{R}}$ and ${}_w^{c3}{\boldsymbol{R}}$ represent the rotation matrix between the coordinate systems of each camera, respectively, and ${}_w^{c1}{\unicode{x1D71E}}$ , ${}_w^{c2}{\unicode{x1D71E}}$ and ${}_w^{c3}{\unicode{x1D71E}}$ denote the translation vectors between the coordinate systems of each camera. The homogeneous coordinates of p in image coordinate system denote ${P_{c1}}{\rm{ }}\left( {{X_{c1}},{Y_{c1}},1} \right)$ , ${P_{c2}}{\rm{ }}\left( {{X_{c2}},{Y_{c2}},1} \right)$ and ${P_{c3}}{\rm{ }}\left( {{X_{c3}},{Y_{c3}},1} \right)$ , respectively.

(25) \begin{align}\left\{ \begin{array}{c}{z_{c1}}{{\boldsymbol{P}}_{c1}}^{\rm{T}} = {z_{c1}} \left[ {\begin{array}{*{20}{c}}{{X_{c1}}}\\ {{Y_{c1}}}\\ 1\end{array}} \right] = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} {{f_{c1}}}& {}0& {}0& {}0\\ 0& {}{{f_{c1}}}& {}0& {}0\\ 0& {}0 {}&1& {}0\end{array}} \right]\left[ {\begin{array}{*{20}{c}}{{x_{c1}}}\\ {{y_{c1}}}\\ {{z_{c1}}}\\ 1\end{array}} \right] = \left[ {{{\boldsymbol{F}}_{c1}}{\rm{ }}0} \right]{{\boldsymbol{P}}_{c1}}^{\rm{T}}\\[23pt] {z_{c2}}{{\boldsymbol{P}}_{c2}}^{\rm{T}} = {z_{c2}}\left[ {\begin{array}{*{20}{c}}{{X_{c2}}}\\ {{Y_{c2}}}\\1\end{array}} \right] = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} {{f_{c2}}}& {}0& {}0 {}&0\\ 0 {}&{{f_{c2}}}& {}0& {}0\\ 0& {}0 {}&1 {}&0\end{array}} \right]\left[ {\begin{array}{*{20}{c}}{{x_{c2}}}\\ {{y_{c2}}}\\ {{z_{c2}}}\\ 1\end{array}} \right] = \left[ {{{\boldsymbol{F}}_{c2}}{\rm{ }}0} \right]{{\boldsymbol{P}}_{c2}}^{\rm{T}}\\[23pt] {z_{c3}}{{\boldsymbol{P}}_{c3}}^{\rm{T}} = {z_{c3}}\left[ {\begin{array}{*{20}{c}}{{X_{c3}}}\\ {{Y_{c3}}}\\1\end{array}} \right] = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} {{f_{c3}}}& {}0& {}0 {}&0\\ 0& {}{{f_{c3}}}& {}0& {}0\\ 0& {}0& {}1& {}0\end{array}} \right]\left[ {\begin{array}{*{20}{c}}{{x_{c3}}}\\ {{y_{c3}}}\\ {{z_{c3}}}\\ 1\end{array}} \right] = \left[ {{{\boldsymbol{F}}_{c3}}{\rm{ }}0} \right]{{\boldsymbol{P}}_{c3}}^{\rm{T}}\end{array} \right. \end{align}

where ${f_{c1}}$ , ${f_{c2}}$ , and ${f_{c3}}$ denote the focal lengths of the camera 1, camera 2 and camera 3, respectively. The homogeneous coordinates of pixels corresponding to its imaging point denote ${p_1}\left( {{u_1},{v_1},1} \right)$ , ${p_2}\left( {{u_2},{v_2},1} \right)$ and ${p_3}\left( {{u_3},{v_3},1} \right)$ , respectively.

(26) \begin{align}\left\{ \begin{array}{c} {\boldsymbol{P}}_{1}^{T} = \left[ \begin{array}{c} u_{1}\\ v_{1}\\ 1 \end{array} \right] = \left[ \begin{array}{*{20}{ccc}} 1/dX_{c1} & \gamma & {}u_{01}\\ 0 & {}1/dY_{c1} & {}v_{01}\\ 0& {}0 {}& 1 \end{array} \right] = \left[ \begin{array}{*{20}{c}} X_{1}\\ Y_{1}\\ 1 \end{array} \right] = {\boldsymbol{A}}_{c1} {\boldsymbol{P}}_{c1}^{\rm{T}}\\[15pt] {\boldsymbol{P}}_{3}^{T} = \left[ \begin{array}{c} u_{3}\\ v_{3}\\ 1 \end{array} \right] = \left[ \begin{array}{*{20}{ccc}} 1/dX_{c3} & \gamma & {}u_{03}\\ 0 & {}1/dY_{c3} & {}v_{03}\\ 0& {}0 {}& 1 \end{array} \right] = \left[ \begin{array}{*{20}{c}} X_{3}\\ Y_{3}\\ 1 \end{array} \right] = {\boldsymbol{A}}_{c3} {\boldsymbol{P}}_{c3}^{\rm{T}}\\[15pt] {\boldsymbol{P}}_{3}^{T} = \left[ \begin{array}{c} u_{3}\\ v_{3}\\ 1 \end{array} \right] = \left[ \begin{array}{*{20}{ccc}} 1/dX_{c3} & \gamma & {}u_{03}\\ 0 & {}1/dY_{c3} & {}v_{03}\\ 0& {}0 {}& 1 \end{array} \right] = \left[ \begin{array}{*{20}{c}} X_{3}\\ Y_{3}\\ 1 \end{array} \right] = {\boldsymbol{A}}_{c3} {\boldsymbol{P}}_{c3}^{\rm{T}} \end{array} \right. \end{align}

where $d{X_{c1}}d{Y_{c1}}$ , $d{X_c}_2d{Y_{c2}}$ , $d{X_c}_3d{Y_{c3}}$ represent the lengths of each row and column of each camera pixel, respectively. $\gamma $ denotes the slope of the corresponding axis. $\left( {{u_{01}},{v_{01}}} \right)$ , $\left( {{u_{02}},{v_{02}}} \right)$ and $\left( {{u_{03}},{v_{03}}} \right)$ denote the origin coordinates of the image coordinate system of each camera on the pixel coordinate system. Integrating Equations (24) to (26), the transformation relationship between point p and its corresponding points ${p_1}$ , ${p_2}$ and ${p_3}$ in each camera pixel coordinates can be obtained as follows.

(27) \begin{align}\left\{ \begin{array}{c}s{{\boldsymbol{p}}_1}^{\rm{T}} = s\left[ {\begin{array}{*{20}{c}}{{u_1}}\\ {{v_1}}\\ 1\end{array}} \right] = \left[ {\begin{array}{*{20}{l}}{{f_{{X_{c1}}}}}&\gamma&{{u_{01}}}\\ 0&{{f_{{Y_{c1}}}}}&d{{v_{01}}}\\ 0&0&1\end{array}} \right] \left[ {\begin{array}{*{20}{c}}{{}_w^{c1}{r_{11}}} {}\quad {{}_w^{c1}{r_{12}}}\quad {}{{}_w^{c1}{r_{13}}} {}\quad{{}_w^{c1}{\tau _x}}\\[3pt] {{}_w^{c1}{r_{21}}} {} \quad{{}_w^{c1}{r_{22}}}\quad {}{{}_w^{c1}{r_{23}}} {}\quad{{}_w^{c1}{\tau _y}}\\[3pt] {{}_w^{c1}{r_{31}}} {} \quad{{}_w^{c1}{r_{32}}} {}\quad{{}_w^{c1}{r_{33}}}\quad {}{{}_w^{c1}{\tau _z}}\end{array}} \right] \left[ {\begin{array}{*{20}{c}}{{x_w}}\\ {{y_w}}\\ {{z_w}}\\ 1\end{array}} \right] = {{\boldsymbol{D}}_{c1}}{}_w^{c1}{\boldsymbol{M}}{\boldsymbol{p}^{\rm{T}}}\\[28pt] s{{\boldsymbol{p}}_2}^{\rm{T}} = s\left[ {\begin{array}{*{20}{c}}{{u_2}}\\ {{v_2}}\\ 1\end{array}} \right] = \left[ {\begin{array}{*{20}{l}}{{f_{{X_{c2}}}}}&\gamma&{{u_{02}}}\\ 0&{{f_{{Y_{c2}}}}}&{{v_{02}}}\\ 0&0&1\end{array}} \right] \left[ {\begin{array}{*{20}{l}}{{}_w^{c2}{r_{11}}}\quad {}{{}_w^{c2}{r_{12}}}\quad {}{{}_w^{c2}{r_{13}}}\quad {}{{}_w^{c2}{\tau _x}}\\[3pt] {{}_w^{c2}{r_{21}}}\quad {}{{}_w^{c2}{r_{22}}}\quad {}{{}_w^{c2}{r_{23}}}\quad {}{{}_w^{c2}{\tau _y}}\\[3pt] {{}_w^{c2}{r_{31}}}\quad {}{{}_w^{c2}{r_{32}}}\quad {}{{}_w^{c2}{r_{33}}}\quad {}{{}_w^{c2}{\tau _z}}\end{array}} \right] \left[ {\begin{array}{*{20}{c}}{{x_w}}\\ {{y_w}}\\ {{z_w}}\\ 1\end{array}} \right] = {{\boldsymbol{D}}_{c2}}{}_w^{c2}{\boldsymbol{M}}{{\boldsymbol{p}}^{\rm{T}}}\\[28pt] s{{\boldsymbol{p}}_3}^{\rm{T}} = s\left[ {\begin{array}{*{20}{c}}{{u_3}}\\ {{v_3}}\\ 1\end{array}} \right] = \left[ {\begin{array}{*{20}{l}}{{f_{{X_{c3}}}}}&\gamma&{{u_{03}}}\\ 0&{{f_{{Y_{c3}}}}}&{{v_{03}}}\\ 0&0&1\end{array}} \right] \left[ {\begin{array}{*{20}{l}}{{}_w^{c3}{r_{11}}}\quad {}{{}_w^{c3}{r_{12}}}\quad {}{{}_w^{c3}{r_{13}}} {}\quad{{}_w^{c3}{\tau _x}}\\[3pt] {{}_w^{c3}{r_{21}}}\quad {}{{}_w^{c3}{r_{22}}}\quad {}{{}_w^{c3}{r_{23}}}\quad {}{{}_w^{c3}{\tau _y}}\\[3pt] {{}_w^{c3}{r_{31}}} {}\quad {{}_w^{c3}{r_{32}}}\quad {}{{}_w^{c3}{r_{33}}}\quad {}{{}_w^{c3}{\tau _z}}\end{array}} \right] \left[ {\begin{array}{*{20}{c}}{{x_w}}\\ {{y_w}}\\ {{z_w}}\\ 1\end{array}} \right] = {{\boldsymbol{D}}_{c3}}{}_w^{c3}{\boldsymbol{M}}{{\boldsymbol{p}}^{\rm{T}}}\end{array} \right. \end{align}

where s denotes nonzero proportionality factor. $f_{X_{c}}= f_{c}/dX_{c} \cdot f_{Y_{c}} = f_{c}/dY_{c} \cdot $ c₁, c₂, and c₃ denote Camera 1, Camera 2, and Camera 3, respectively. ${{\boldsymbol{D}}_{c1}}$ , ${{\boldsymbol{D}}_{c2}}$ , ${{\boldsymbol{D}}_{c3}}$ denote intrinsic parameter matrices of camera 1, camera 2 and camera 3. ${}_w^{c1}{\boldsymbol{M}}$ , ${}_w^{c2}{\boldsymbol{M}}$ , ${}_w^{c3}{\boldsymbol{M}}$ denote extrinsic parameter matrices of camera 1, camera 2 and camera 3.

4.2 Multi-cameras vision parameter model

The visual parameters between the camera coordinate systems are used to represent the pose relationship between the cameras, as shown in Fig. 6.

Figure 6.

Visual parameters between camera coordinate systems.

The transformation matrices between each camera are constructed by utilising point p in conjunction with the coordinate transformations between points pc1, pc2 and pc3.

(28) \begin{align}\left\{ \begin{array}{c}{\boldsymbol{p}}_{c2}^{\mathop{\rm T}\nolimits} = \left[ {\begin{array}{*{20}{c}}{{}_{c1}^{c2}{r_{11}}} \quad {}{{}_{c1}^{c2}{r_{12}}} \quad {}{{}_{c1}^{c2}{r_{13}}} {}\quad {{}_{c1}^{c2}{\tau _x}}\\[3pt] {{}_{c1}^{c2}{r_{21}}} \quad {}{{}_{c1}^{c2}{r_{22}}} {}\quad{{}_{c1}^{c2}{r_{23}}} {}\quad{{}_{c1}^{c2}{\tau _y}}\\[3pt] {{}_{c1}^{c2}{r_{31}}} \quad {}{{}_{c1}^{c2}{r_{32}}} {}\quad{{}_{c1}^{c2}{r_{33}}}\quad {}{{}_{c1}^{c2}{\tau _z}}\end{array}} \right] \left[ {\begin{array}{*{20}{c}}{{x_{c1}}}\\ {{y_{c1}}}\\ {{z_{c1}}}\\ 1\end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{}_{c1}^{c2}{R}} {}\quad{{}_{c1}^{c2}{\unicode{x1D71E}}} \end{array}} \right] {\boldsymbol{p}}_{c1}^{\mathop{\rm T}\nolimits} = {}_{c1}^{c2}{M}{\boldsymbol{p}}_{c1}^{\mathop{\rm T}\nolimits} \\[28pt] {\boldsymbol{p}}_{c3}^{\mathop{\rm T}\nolimits} = \left[ {\begin{array}{*{20}{c}} {{}_{c2}^{c3}{r_{11}}} \quad {}{{}_{c2}^{c3}{r_{12}}} \quad {}{{}_{c2}^{c3}{r_{13}}}\quad {}{{}_{c2}^{c3}{\tau _x}}\\[3pt] {{}_{c2}^{c3}{r_{21}}}\quad {}{{}_{c2}^{c3}{r_{22}}}\quad {}{{}_{c2}^{c3}{r_{23}}}\quad {}{{}_{c2}^{c3}{\tau _y}}\\[3pt] {{}_{c2}^{c3}{r_{31}}}\quad {}{{}_{c2}^{c3}{r_{32}}}\quad {}{{}_{c2}^{c3}{r_{33}}}\quad {}{{}_{c2}^{c3}{\tau _z}}\end{array}} \right] \left[ {\begin{array}{*{20}{c}}{{x_{c2}}}\\ {{y_{c2}}}\\ {{z_{c2}}}\\ 1\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}{{}_{c2}^{c3}{R}}\quad {}{{}_{c2}^{c3}{\unicode{x1D71E} }}\end{array}} \right] {\boldsymbol{p}}_{c1}^{\mathop{\rm T}\nolimits} = {}_{c2}^{c3}{M}{\boldsymbol{p}}_{c2}^{\mathop{\rm T}\nolimits} \\[28pt] {\boldsymbol{p}}_{c1}^{\mathop{\rm T}\nolimits} = \left[ {\begin{array}{*{20}{c}} {{}_{c3}^{c1}{r_{11}}}\quad {}{{}_{cr}^{c1}{r_{12}}} {}\quad{{}_{cr}^{c1}{r_{13}}} {}\quad{{}_{cr}^{c1}{\tau _x}}\\[3pt] {{}_{c3}^{c1}{r_{21}}}\quad {}{{}_{c3}^{c1}{r_{22}}} {}\quad{{}_{c3}^{c1}{r_{23}}} {}\quad{{}_{c3}^{c1}{\tau _y}}\\[3pt] {{}_{c3}^{c1}{r_{31}}}\quad {}{{}_{c3}^{c1}{r_{32}}} {}\quad{{}_{c3}^{c1}{r_{33}}} {}\quad{{}_{c3}^{c1}{\tau _z}}\end{array}} \right]\left[ {\begin{array}{*{20}{c}}{{x_{c3}}}\\ {{y_{c3}}}\\ {{z_{c3}}}\\ 1\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}{{}_{c3}^{c1}{R}} \quad {}{{}_{c3}^{c1}{\unicode{x1D71E} }}\end{array}} \right]{\boldsymbol{p}}_{c1}^{\mathop{\rm T}\nolimits} = {}_{c3}^{c1}{M}{\boldsymbol{p}}_{c3}^{\mathop{\rm T}\nolimits} \end{array} \right. \end{align}

where ${}_{c1}^{c2}{\boldsymbol{R}}$ , ${}_{c2}^{c3}{\boldsymbol{R}}$ , ${}_{c3}^{c1}{\boldsymbol{R}}$ represent rotation matrices between each camera, respectively. ${}_{c1}^{c2}{\unicode{x1D71E} }$ , ${}_{c2}^{c3}{\unicode{x1D71E} }$ , ${}_{c3}^{c1}{\unicode{x1D71E} }$ represent translation vectors between each camera, respectively. ${}_{c1}^{c2}{\boldsymbol{M}}$ , ${}_{c2}^{c3}{\boldsymbol{M}}$ , ${}_{c3}^{c1}{\boldsymbol{M}}$ represent homogeneous transformation matrices between each camera, respectively.

5.1 Intrinsic and extrinsic parameters solution

Based on the calibration board, the world coordinate system ${o_w} - {x_w}{y_w}{z_w}$ is established. The calibration board is on the ${z_w} = 0$ , as shown in Fig. 7.

Figure 7.

Multi-cameras calibration.

In Fig. 7, p represents the feature points on the calibration board. The calibration board can be placed arbitrarily within the field of view of camera 1. Therefore, the establishment of the world coordinate system is not unique, even can be arbitrarily established. In order to facilitate the solution, the world coordinate system can be set in the opposite direction to camera 1. Based on it, Equation (27) is rewritten as:

(29) \begin{align}\left\{ \begin{array}{c}s{{\boldsymbol{p}}_1}^{\rm{T}} = {{\boldsymbol{D}}_{c1}} \left[ {\begin{array}{*{20}{c}} {{}_w^{c1}{{\boldsymbol{r}}_1}} \quad {}{{}_w^{c1}{{\boldsymbol{r}}_2}}\quad {}{{}_w^{c1}{{\boldsymbol{r}}_3}}\quad {}{{}_w^{c1}{\bf{\tau }}}\end{array}}\! \right] \left[ {\begin{array}{*{20}{c}}{{x_w}}\\ {{y_w}}\\ 0\\ 1\end{array}} \right] = {{\boldsymbol{D}}_{c1}} \left[ {\begin{array}{*{20}{c}}{{}_w^{c1}{{\boldsymbol{r}}_1}}\quad {}{{}_w^{c1}{{\boldsymbol{r}}_2}}\quad {}{{}_w^{c1}{\bf{\tau }}}\end{array}} \!\right] \left[ {\begin{array}{*{20}{c}}{{x_w}}\\ {{y_w}}\\ 1\end{array}} \right]\\[28pt] s{{\boldsymbol{p}}_2}^{\rm{T}} = {{\boldsymbol{D}}_{c2}} \left[ {\begin{array}{*{20}{c}}{{}_w^{c2}{{\boldsymbol{r}}_1}}\quad {}{{}_w^{c2}{{\boldsymbol{r}}_2}}\quad {}{{}_w^{c2}{{\boldsymbol{r}}_3}}\quad {}{{}_w^{c2}{\bf{\tau }}}\end{array}}\! \right] \left[ {\begin{array}{*{20}{c}}{{x_w}}\\ {{y_w}}\\ 0\\ 1\end{array}} \right] = {{\boldsymbol{D}}_{c2}} \left[ {\begin{array}{*{20}{c}}{{}_w^{c2}{{\boldsymbol{r}}_1}} {}\quad{{}_w^{c2}{{\boldsymbol{r}}_2}}\quad {}{{}_w^{c2}{\bf{\tau }}}\end{array}}\! \right] \left[ {\begin{array}{*{20}{c}}{{x_w}}\\ {{y_w}}\\ 1\end{array}} \right]\\[28pt] s{{\boldsymbol{p}}_3}^{\rm{T}} = {{\boldsymbol{D}}_{c3}} \left[ {\begin{array}{*{20}{c}}{{}_w^{c3}{{\boldsymbol{r}}_1}} {}\quad{{}_w^{c3}{{\boldsymbol{r}}_2}} {}\quad{{}_w^{c3}{{\boldsymbol{r}}_3}}\quad {}{{}_w^{c3}{\bf{\tau }}}\end{array}}\! \right] \left[ {\begin{array}{*{20}{c}}{{x_w}}\\ {{y_w}}\\ 0\\ 1\end{array}} \right] = {{\boldsymbol{D}}_{c3}} \left[ {\begin{array}{*{20}{c}}{{}_w^{c3}{{\boldsymbol{r}}_1}}\quad {}{{}_w^{c3}{{\boldsymbol{r}}_2}}\quad {}{{}_w^{c3}{\bf{\tau }}}\end{array}} \!\right] \left[ {\begin{array}{*{20}{c}}{{x_w}}\\ {{y_w}}\\ 1\end{array}} \right]\end{array} \right. \end{align}

Then there exists a homography matrix.

(30) \begin{align}{\boldsymbol{H}} = \left[ {\begin{array}{*{20}{c}}{{{\boldsymbol{h}}_1}} {}\quad{{{\boldsymbol{h}}_2}} \quad {}{{{\boldsymbol{h}}_3}}\end{array}} \right] = \lambda {{\boldsymbol{D}}_c}\left[ {\begin{array}{*{20}{c}}{{}_w^c{{\boldsymbol{r}}_1}} {} \quad{{}_w^c{{\boldsymbol{r}}_2}} {}\quad{{}_w^c{\bf{\tau }}}\end{array}} \right] \end{align}

where $\lambda $ is an arbitrary scalar c represent camera 1, camera 2 and camera 3. The homography matrix can be represented as the mapping relationship between feature point p on the calibration board and its corresponding pixel coordinates on camera 1, camera 2 and camera 3. According to the orthogonality of ${}_w^c{\boldsymbol{R}}$ , ${}_w^c{{\boldsymbol{r}}_1}^{\rm{T}}{}_w^c{{\boldsymbol{r}}_2} = 0$ and ${}_w^c{{\boldsymbol{r}}_1}^{\rm{T}}{}_w^c{{\boldsymbol{r}}_1} = {}_w^c{{\boldsymbol{r}}_2}^{\rm{T}}{}_w^c{{\boldsymbol{r}}_2}$ are obtained. Then Equation (31) is obtained.

(31) \begin{align}{\boldsymbol{B}} = {{\boldsymbol{D}}_c}^{ - {\rm{T}}}{{\boldsymbol{D}}_c}^{ - 1} \!=\! \left[ {\begin{array}{c@{\kern5pt}c@{\kern5pt}c} b_{11} & b_{12} & b_{13}\\ b_{21} & b_{22} & b_{23}\\ b_{31} & b_{32} & b_{33}\end{array}} \right] \left. {\! =\! \left[ {\begin{array}{*{20}{c}}{\dfrac{1}{{f_{{X_c}}^2}}} {}&{ - \dfrac{\gamma }{{f_{{X_c}}^2{\,f_{{Y_c}}}}}} & {}{\dfrac{{{v_0}\gamma - {u_0}{\,f_{{Y_c}}}}}{{f_{{X_c}}^2{\,f_{{Y_c}}}}}}\\ { - \dfrac{\gamma }{{f_{{X_c}}^2{\,f_{{Y_c}}}}}}& {}{\dfrac{{{\gamma ^2}}}{{\,f_{{X_c}}^2f_{{Y_c}}^2}} + \dfrac{1}{{f_{{Y_c}}^2}}} {} &{ - \dfrac{{\gamma \left( {{v_0}\gamma - {u_0}{\,f_{{Y_c}}}} \right)}}{{f_{{X_c}}^2{\,f_{{Y_c}}}}} - \dfrac{{{v_0}}}{{f_{{Y_c}}^2}}}\\ {\dfrac{{{v_0}\gamma - {u_0}{\,f_{{Y_c}}}}}{{f_{{X_c}}^2{\,f_{{Y_c}}}}}} {} &{ - \dfrac{{\gamma \left( {{v_0}\gamma - {u_0}{\,f_{{Y_c}}}} \right)}}{{f_{{X_c}}^2{\,f_{{Y_c}}}}} - \dfrac{{{v_0}}}{{f_{{Y_c}}^2}}}& {}{\dfrac{{{{\left( {{v_0}\gamma - {u_0}{\,f_{{Y_c}}}} \right)}^2}}}{{\,f_{{X_c}}^2\,f_{{Y_c}}^2}} + \dfrac{{v_0^2}}{{f_{{Y_c}}^2}} + 1}\end{array}} \right.} \right] \end{align}

Let the i-th column vector of matrix ${\boldsymbol{H}}$ be ${{\boldsymbol{h}}_i} = {\left[ {\begin{array}{*{20}{c}}{{h_{1i}}} {}{{h_{2i}}} {}{{h_{3i}}}\end{array}} \right]^{\rm{T}}}$ . j-th column is the same as above. It can be seen that the Equation (31) is a symmetric matrix. Then a six-dimensional vector ${\boldsymbol{b}} = {\left[ {\begin{array}{*{20}{c}}{{b_{11}}} \quad {}{{b_{12}}} {} \quad{{b_{13}}} {} \quad{{b_{22}}}\quad {}{{b_{23}}}\quad {}{{b_{33}}}\end{array}} \right]^{\rm{T}}}$ is defined. There exists a vector ${{\boldsymbol{E}}_{ij}}$ such that ${\boldsymbol{h}}_i^{\rm{T}}{\boldsymbol{B}}{{\boldsymbol{h}}_j} = {\boldsymbol{E}}_{ij}^{\rm{T}}{\boldsymbol{b}}$ and ${{\boldsymbol{E}}_{ij}} = \left[ {\begin{array}{*{20}{c}}{{h_{1i}}{h_{1j}}} \quad {}{{h_{1i}}{h_{2j}} + {h_{2i}}{h_{1j}}}\quad {}{{h_{2i}}{h_{3j}}} {}\quad{{h_{3i}}{h_{1j}} + {h_{1i}}{h_{3j}}} \quad {}{{h_{3i}}{h_{2j}} + {h_{2i}}{h_{3j}}} {} \quad{{h_{3i}}{h_{3j}}}\end{array}} \right]$ . Therefore, it can be seen that any calibration image used can yield two equations as shown below.

(32) \begin{align}\left[ {\begin{array}{*{20}{c}}{{\boldsymbol{E}}_{12}^{\rm{T}}}\\{{{\left( {{{\boldsymbol{E}}_{11}} - {{\boldsymbol{E}}_{22}}} \right)}^{\rm{T}}}}\end{array}} \right]{\boldsymbol{b}} = 0 \end{align}

According to n calibration images taken from different directions, 2n equations can be obtained, which can be expressed as ${\boldsymbol{Eb}} = 0$ . The eigenvector corresponding to the smallest eigenvalue of ${{\boldsymbol{E}}^{\rm{T}}}{\boldsymbol{E}}$ is the least squares solution, which allows for the determination of the six-dimensional vector b . Consequently, the camera intrinsic parameters can also be solved.

(33) \begin{align} \left\{ \begin{array}{l} {v_0} = \left(b_{12}b_{13} - b_{1}b_{23}\right) / (b_{11}b_{22} - b_{12}^{2})\\[4pt] \lambda = b_{33} - \left[b_{13}^2 + v_{0}(b_{12}b_{13} - b_{11}b_{23})\right]/b_{11}\\[4pt] f_{X_{c}} = \sqrt{\lambda/b_{11}}\\[4pt] f_{Y_{c}} = \sqrt{\lambda/b_{11}(b_{11}b_{22}-b_{12}^{2})}\\[4pt] \lambda = -b_{12}f^{2}_{X_{c}}f_{Y_{c}/\lambda}\\[4pt] u_{0} = \lambda v_{0}/f_{X_{c}}-b_{13}f^{2}_{X_{c}}/\lambda \end{array} \right. \end{align}

Therefore, extrinsic parameters can be obtained.

(34) \begin{align}\left\{ \begin{array}{l}{}_w^c{{\boldsymbol{r}}_1} = \lambda {\boldsymbol{D}}_c^{ - 1}{{\boldsymbol{h}}_1}\\[3pt] {}_w^c{{\boldsymbol{r}}_2} = \lambda {\boldsymbol{D}}_c^{ - 1}{{\boldsymbol{h}}_2}\\[3pt] {}_w^c{{\boldsymbol{r}}_3} = {{\boldsymbol{r}}_1} \times {{\boldsymbol{r}}_2}\\[3pt] {}_w^c{\bf{\tau }} = \lambda {\boldsymbol{D}}_c^{ - 1}{{\boldsymbol{h}}_3}\end{array} \right.\quad \end{align}

In Equation (34), $\lambda = {{1{\rm{ }}} \mathord{\left/ {{\rm{ }}\left\| {{\boldsymbol{D}}_c^{ - 1}{{\boldsymbol{h}}_1}} \right\|} \right.} = {{1{\rm{ }}} \mathord{\left/ {{\rm{ }}\left\| {{\boldsymbol{D}}_c^{ - 1}{{\boldsymbol{h}}_2}} \right\|} \right. }} }$ . Due to the influence of noise, the optimal rotation matrix needs to be obtained by singular value decomposition and other methods.

5.2 Solution of visual parameter between cameras

The pose relationship between the camera coordinate systems can be obtained by Equation (28). Combined with Equation (25), the mapping relationship between the two imaging points of the feature point p in the image coordinate systems of camera 1 and camera 2 can be solved.

(35) \begin{align}\begin{array}{c}{z_{c2}}\left[ {\begin{array}{*{20}{c}}{{X_{c2}}}\\ {{Y_{c2}}}\\1\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}{{f_{c2}}}\quad {}0\quad {}0\\ 0\quad {}{{f_{c2}}}\quad {}0\\ 0\quad {}0\quad {}1\end{array}} \right]\left[ {\begin{array}{*{20}{c}}{{}_{c1}^{c2}{\boldsymbol{R}}}\quad {}{{}_{c1}^{c2}{\unicode{x1D71E}}}\end{array}} \right] \left[ {\begin{array}{*{20}{c}}{{x_{c1}}}\\ {{y_{c1}}}\\ {{z_{c1}}}\\ 1\end{array}} \right]\\[23pt] = \left[ {\begin{array}{{c@{\quad}c@{\quad}c@{\quad}c}}{{f_{c2}}{}_{c1}^{c2}{r_{11}}} & {{f_{c2}}{}_{c1}^{c2}{r_{12}}} & {{f_{c2}}{}_{c1}^{c2}{r_{13}}} & {{f_{c2}}{}_{c1}^{c2}{\tau _x}}\\[4pt] {{f_{c2}}{}_{c1}^{c2}{r_{21}}} & {{f_{c2}}{}_{c1}^{c2}{r_{22}}} & {{f_{c2}}{}_{c1}^{c2}{r_{23}}} & {{f_{c2}}{}_{c1}^{c2}{\tau _y}}\\[4pt] {{}_{c1}^{c2}{r_{31}}} & {{}_{c1}^{c2}{r_{32}}} & {{}_{c1}^{c2}{r_{33}}} & {{}_{c1}^{c2}{\tau _z}}\end{array}} \right] \left[ {\begin{array}{*{20}{c}}{\dfrac{{{X_{c1}}{z_{c1}}}}{{{f_{c1}}}}}\\[8pt] {\dfrac{{{Y_{c1}}{z_{c1}}}}{{{f_{c1}}}}}\\{{z_{c1}}}\\[8pt] 1\end{array}} \right]\end{array} \end{align}

The mapping relationship between any other two camera imaging points can also be solved by the method described above. According to Equation (35), Equation (36) is obtained to represent the relationship between the three-dimensional coordinates ${p_w}\left( {{x_w},{y_w},{z_w}} \right)$ of the feature point p in the world coordinate system and other parameters.

(36) \begin{align}\left\{ \begin{array}{l}{x_w} = \dfrac{{{X_{c1}}{z_w}}}{{{f_{c1}}}}\\[10pt] {y_w} = \dfrac{{{Y_{c1}}{z_w}}}{{{f_{c1}}}}\\ {z_w} = \dfrac{{{f_{c1}}\left( {{f_{c2}}{}_{c1}^{c2}{\tau _x} - {X_{c2}}{}_{c1}^{c2}{\tau _z}} \right)}} {{{X_{c2}}\left( {{}_{c1}^{c2}{r_{31}}{X_{c1}} + {}_{c1}^{c2}{r_{32}}{Y_{c1}} + {f_{c1}}{}_{c1}^{c2}{r_{33}}} \right) - {f_{c2}}\left( {{}_{c1}^{c2}{r_{11}}{X_{c1}} + {}_{c1}^{c2}{r_{12}}{Y_{c1}} + {f_{c1}}{}_{c1}^{c2}{r_{13}}} \right)}}\end{array} \right. \end{align}

After obtaining the coordinates of multiple feature points on the calibration board taken in different directions, the visual parameters between each camera can be solved by incorporating the intrinsic parameters.

5.3 Distortion correction

The aforementioned process of solving for intrinsic and extrinsic parameters does not take into account the distortion of camera lens. In practical applications, the camera lens will have certain geometric distortion, including radial distortion and tangential distortion. In the calibration process, what is known is the image that has undergone distortion. In order to obtain the image without distortion, it is necessary to derive the mapping relationship through the distortion correction model.

Due to the distortions introduced by camera lens, there is a certain deviation between the image point of point p with distortion and its undistorted counterpart. This deviation also affects the accuracy of deviation detection of fork-ear. Therefore, to enhance the detection accuracy of the camera, it is necessary to consider the radial distortions and tangential distortions of the lens.

The offset coordinate of the actual imaging point p from $p\left( {u,v} \right)$ to $p'\left( {{u_d},{v_d}} \right)$ is as follows.

(37) \begin{align}\left\{ \begin{array}{l}{u_d} = u + \Delta u\\{v_d} = v + \Delta v\end{array} \right. \end{align}

1. (1) Radial distortion model

The radial distortion correction model is given by

(38) \begin{align}\left\{ {\begin{array}{*{20}{c}}{u' = u(1 + {k_1}{\mu ^2} + {k_2}{\mu ^4})}\\[4pt] {v' = v(1 + {k_1}{\mu ^2} + {k_2}{\mu ^4})}\end{array}} \right. \end{align}

where $\left( {u',v'} \right)$ denotes the pixel coordinates of the distorted image point under the distortion condition. $\left( {u,v} \right)$ denotes the pixel coordinates of the undistorted image point in the pinhole camera model, ${\mu ^2} = {u^2} + {v^2}$ . ${k_1}$ and ${k_2}$ are the radial distortion coefficients.

1. (2) Tangential distortion model

The tangential distortion correction model is given by

(39) \begin{align}\left\{ {\begin{array}{*{20}{c}}{u' = u + 2{q_1}uv + {q_2}({\mu ^2} + 2{u^2})}\\[4pt] {v' = v + {q_1}({\mu ^2} + 2{v^2}) + 2{q_2}uv}\end{array}} \right. \end{align}

where ${q_1}$ and ${q_2}$ are tangential distortion. By combining Equations (37) and (38), the distortion correction model is obtained.

(40) \begin{align}\left\{ {\begin{array}{*{20}{c}}{u' = u + {k_1}u{r^2} + {k_2}u{r^4} + 2{q_1}u\nu + {q_2}(2{u^2} + {\mu ^2})}\\[4pt] {v' = v + {k_1}v{r^2} + {k_2}v{r^4} + {q_1}({\mu ^2} + 2{u^2}) + 2{q_2}u\nu }\end{array}} \right. \end{align}

The lens distortion model is a nonlinear mathematical model, which cannot be calculated by the direct method. Therefore, the camera calibration parameters of the pinhole model are used as the initial value. Equations (37) and (40) are substituted into Equation (27) to construct the reprojection equation. The lens distortion parameters are optimised by minimising the reprojection error as the objective.

(41) \begin{align}{p_{ij}} = \mathop {\mathop \sum }\limits_{n}^{i = 1} \mathop {\mathop \sum }\limits_{m}^{j = 1} {m_{ij}} - P{(K,D,{R_i},{T_i},{M_j})^2} \end{align}

where D represents the lens distortion parameters. n is the total number of target images captured; m is the number of feature points. ${p_{ij}}$ denotes the image coordinates of the j-th feature point on the i-th image after outlier images are excluded; P is the reprojection equation. ${R_i}$ and ${T_i}$ represent the rotation matrix and translation vector from the camera coordinates to the world coordinate system of the i-th image after outlier images are excluded. ${M_j}$ represents the world coordinates of the j-th circle centre feature point in the image. Ultimately, the Levenberg-Marquardt optimisation algorithm is employed to iteratively solve Equation (41), thereby determining the lens distortion parameters.