2.1. Camera coordinate frame

As shown in Fig. 2, (X _c , Y _c , Z _c ) and (X _w , Y _w , Z _w ) represent the camera and world coordinate frames, respectively. For a point located at the center of the ArUco marker, its camera and world coordinate values are (X _c0 , Y _c0 , Z _c0 ) and (X _w0 , Y _w0 , Z _w0 ), respectively. The marker has a square shape and four corner points, which are denoted by (X _c1 , Y _c1 , Z _c1 ), (X _c2 , Y _c2 , Z _c2 ), (X _c3 , Y _c3 , Z _c3 ), and (X _c4 , Y _c4 , Z _c4 ) in the camera coordinate frame. In our lab experiments, a process (detailed in Section 3) is developed to align the cameras almost parallel to the floorboard. The ArUco marker is placed flat on the floorboard and has a small size; therefore, it is valid to assume $Z_{c1}=Z_{c2}=Z_{c3}=Z_{c4}$ .

Figure 2. Camera and world coordinate frames.

Figure 3. Normalized image coordinates and camera coordinates of the ArUco marker.

The length L of the four edges of the ArUco marker in the camera coordinate is the same and can be expressed as

(1) \begin{align} L&=\sqrt{(X_{cj}-X_{ci})^{2}+(Y_{cj}-Y_{ci})^{2}+(Z_{cj}-Z_{ci})^{2}}\;\approx\; \sqrt{(X_{cj}-X_{ci})^{2}+(Y_{cj}-Y_{ci})^{2}} \nonumber \\[5pt]&= \frac{1}{4} \sum_{i=1}^{4} \sqrt{(X_{cj}-X_{ci})^{2}+(Y_{cj}-Y_{ci})^{2}} \end{align}

where i = 1, 2, 3, and 4, j = mod(i, 4) + 1, and ‘mod’ operation denotes the remainder after division. Usually, the values of the corners (X _ci , Y _ci , Z _ci ) in the camera coordinate frame are unknown. However, we can find them from their normalized coordinates (u, ν, 1) using OpenCV undistortPoints() function, where u=X _c /Z _c and v=Y _c /Z _c . This function corrects lens distortion and normalizes the coordinates of detected points. According to the triangle similarity Theorems (Fig. 3), Eq. (2) can also be expressed as

(2) \begin{equation} L=\frac{Z_{c0}}{4}\sum\nolimits_{i=1}^{4}\sqrt{\left(u_{j}-u_{i}\right)^{2}+\left(\nu _{j}-\nu _{i}\right)^{2}} \end{equation}

where $Z_{c0}$ is the value of $Z_{c}$ at the marker center, and $Z_{c0}=Z_{c1}=Z_{c2}=Z_{c3}=Z_{c4}$ because the marker is small, flat, and almost parallel to the camera. In this work, the marker length L is measured manually, which allows us to calculate $Z_{c0}$ of the marker by

(3) \begin{equation} Z_{c0}=\frac{4L}{\sum _{i=1}^{4}\sqrt{\left(u_{j}-u_{i}\right)^{2}+\left(\nu _{j}-\nu _{i}\right)^{2}}} \end{equation}

Figure 4. The neural network model used in this study.

As the marker is a square, the x and y coordinates of its center, that is, X _c0 and Y _c0, can be written as

(4) \begin{equation} X_{c0}=\displaystyle\sum_{i=1}^{4}X_{ci}=\frac{Z_{c0}}{4}\sum_{i=1}^{4}u_{i}.\end{equation}

(5) \begin{equation} Y_{c0}=\displaystyle\sum_{i=1}^{4}Y_{ci}=\frac{Z_{c0}}{4}\sum_{i=1}^{4}v_{i}.\end{equation}

Thus, the camera coordinate values of the marker’s center, that is, ( $X_{c0}, Y_{c0}, Z_{c0}$ ) can be completely determined by Eqs. (3)-(5).

2.2. World coordinate frame

The next step is to transform the ArUco marker from the camera coordinate frame to the world coordinate frame for localization in the real environment. To establish the transformation relationship F between them, that is, (X _w , Y _w , Z _w ) = F (X _c , Y _c , Z _c ), the true value of the ArUco marker in the world coordinate frame is needed and can be manually measured by treating one location in the real environment as the origin. Note that the location of the ArUco marker in the camera coordinate is now known following the procedure in Section 2.1. The ArUco marker is placed at multiple locations, and the measurement is repeated accordingly, which yields a dataset containing many pairs of (X _c , Y _c , Z _c ) and (X _w , Y _w , Z _w ), each corresponding to one marker location. The dataset is then split into two groups, respectively, for training and testing of model F .

The transformation relationship F can be identified by various data training/learning approaches, such as rigid transformation, polynomial regression, Kriging interpolation, machine learning, and others, which are described in detail below.

2.2.1. Rigid transformation

Rigid transformation is the most widely used approach to estimate mapping between the camera and the world coordinates. It is a kind of transformation that does not change the Euclidean distance of every pair of points, such as translations and rotations (that are considered in this work). It is mathematically expressed as

(6) \begin{equation}\left[\begin{array}{c} X_{w}\\[5pt] Y_{w}\\[5pt] Z_{w}\\[5pt] 1 \end{array}\right]=\left[\begin{array}{c@{\quad}c} {\boldsymbol{R}}_{3\times 3} & \boldsymbol{\Gamma }_{3\times 1}\\[5pt] \textbf{0}_{1\times 3} & 1 \end{array}\right]\left[\begin{array}{c} X_{c}\\[5pt] Y_{c}\\[5pt] Z_{c}\\[5pt] 1 \end{array}\right]\end{equation}

where R and $\boldsymbol{\Gamma }$ are the rotation matrix and translation matrix, respectively. Given the dataset in the camera and world coordinate frames above, R and $\boldsymbol{\Gamma }$ can be computed by singular value decomposition (SVD) [Reference Arun, Huang and Blostein24, Reference Kurobe, Sekikawa, Ishikawa and Corsnet25]. The data pairs in the camera frame ${\boldsymbol{X}}$ and in the world frame ${\boldsymbol{Y}}$ are

(7) \begin{equation}{\boldsymbol{X}}=\left(\begin{array}{c} {\boldsymbol{x}}^{(1)}\\[5pt] \vdots \\[5pt] {\boldsymbol{x}}^{(n)} \end{array}\right),{\boldsymbol{Y}}=\left(\begin{array}{c} {\boldsymbol{y}}^{(1)}\\[5pt] \vdots \\[5pt] {\boldsymbol{y}}^{(n)} \end{array}\right)\end{equation}

where $\textbf{X}\in \boldsymbol{\mathfrak{R}}^{n\times 3}$ is the matrix composed of n observations for the three measured input quantities ( $X_{c},$ $Y_{c}$ , and $Z_{c}$ ), that is, each entry ${\boldsymbol{x}}^{(i)}=(X_{c}^{(i)},Y_{c}^{(i)}, Z_{c}^{(i)})$ is the ith observation and $0\leq i\leq n$ ; $\textbf{Y}\in \boldsymbol{\mathfrak{R}}^{n\times 3}$ is the matrix containing n measurements for the three output quantities ( $X_{w},$ $Y_{w}$ , or $Z_{w}$ ), and each entry ${\boldsymbol{y}}^{(i)}=(X_{w}^{(i)},Y_{w}^{(i)}, Z_{w}^{(i)})$ . The centroids of X and Y are calculated by

(8) \begin{equation} \overline{{\boldsymbol{X}}}=\frac{1}{n}\sum\nolimits_{i=1}^{n}{\boldsymbol{x}}^{\left(i\right)},\; \textrm{and}\; \overline{{\boldsymbol{Y}}}=\frac{1}{n}\sum\nolimits_{i=1}^{n}{\boldsymbol{y}}^{\left(i\right)} \end{equation}

Note that $\overline{{\boldsymbol{X}}}\in \boldsymbol{\mathfrak{R}}^{1\times 3}$ and $\overline{{\boldsymbol{Y}}}\in \boldsymbol{\mathfrak{R}}^{1\times 3}$ are row vectors in 3D. Then, the covariance matrix H [Reference Arun, Huang and Blostein24, Reference Kurobe, Sekikawa, Ishikawa and Corsnet25] is defined as

(9) \begin{equation} {\boldsymbol{H}}=({\boldsymbol{X}}-\overline{{\boldsymbol{X}}})^{T}({\boldsymbol{Y}}-\overline{{\boldsymbol{Y}}}) \end{equation}

$({\boldsymbol{X}}-\overline{{\boldsymbol{X}}})$ is an operation that each row of ${\boldsymbol{X}}$ is subtracted by $\overline{{\boldsymbol{X}}}$ . Then, H is decomposed by SVD to produce U and V :

(10) \begin{equation} [{\boldsymbol{U}},{\boldsymbol{S}},{\boldsymbol{V}}]=SVD({\boldsymbol{H}}) \end{equation}

Thus, the R and $\boldsymbol{\Gamma }$ can be obtained

(11) \begin{equation} {\boldsymbol{R}}={\boldsymbol{V}}{\boldsymbol{U}}^{T} \end{equation}

(12) \begin{equation} \boldsymbol{\Gamma }=\overline{{\boldsymbol{Y}}}^{T}-{\boldsymbol{R}}\overline{{\boldsymbol{X}}}^{T} \end{equation}

However, even with camera calibration, the issue of image distortion is still present, leading to a poor estimation of the marker location in the camera coordinate frame. Therefore, rigid transformation solely may be insufficient to yield an accurate estimation of the world coordinate values in the testing stage. Thus, other data-driven modeling methods that can incorporate additional nonlinearity into the mapping relationship F are also examined in this work.

2.2.2. Polynomial regression

Polynomial regression is a widely used approach to build the mapping relationship F between inputs and output responses, in which polynomial terms of inputs are used as regressors. In this work, X _c , Y _c , and Z _c are inputs, and X _w , Y _w , and Z _w are outputs, and F can be determined by

(13) \begin{equation} \textbf{Y}=\boldsymbol{\xi }{\boldsymbol{F}}+{\boldsymbol{v}} \end{equation}

where $\boldsymbol{\xi }\in \boldsymbol{\mathfrak{R}}^{n\times {n_{p}}}$ is a matrix composed of n _p regressors at the n observations. In this work, our regressors include the constant, linear (i.e., $X_{c}$ , $Y_{c}$ , and $Z_{c}$ ), and the second-order nonlinear terms ( $X_{c}Y_{c}, X_{c}Z_{c}, Y_{c}Z_{c}, {X_{c}}^{2}, {Y_{c}}^{2}$ and ${Z_{c}}^{2}$ ) of the input variables, and hence $n_{p}=10.$ F $\in \boldsymbol{\mathfrak{R}}^{{n_{p}}\times 3}$ is a matrix of model coefficients to be estimated, and ${\boldsymbol{v}}\in \boldsymbol{\mathfrak{R}}^{n\times 3}$ is the matrix of the measurement errors for all the three outputs. The ordinary least-squares solution [Reference Sobester, Forrester and Keane26] to F is the best linear unbiased estimator (BLUE) and can be obtained by minimizing a cost function summing over all the squared residuals at each data point, yielding

(14) \begin{equation} \hat{{\boldsymbol{F}}}=(\boldsymbol{\xi }^{T}\boldsymbol{\xi })^{-1}\boldsymbol{\xi }^{T}\textbf{Y} \end{equation}

Once $\hat{{\boldsymbol{F}}}$ is estimated in the offline training stage, it can be used for online real-time estimation of the marker center in the world coordinate frame, that is, X _w0 , Y _w0 , and Z _w0.

2.2.4. Kriging interpolation model

Kriging, proposed by Krige and Sacks, is a data-driven interpolation technique to predict the response surface [Reference Yang, Hong, Wang and Wang28]. It was originally used in geostatistics and gradually applied to various engineering fields. In this work, the Kriging interpolation model is first developed to capture the mapping relationship F between the camera coordinate value, ${\boldsymbol{x}}=(X_{c}, Y_{c},\text{ and }Z_{c})$ and the world coordinate value ${\boldsymbol{y}}=(X_{w}, Y_{w},\text{ or }Z_{w})$ as presented above. The model has two components: a regression model $f({\boldsymbol{x}})$ to estimate the global trend of the data landscape and a Gaussian process model $Z({\boldsymbol{x}})$ with zero mean and variance $\sigma ^{2}$ to capture the difference between the trend function and the true response surface. Therefore, the Kriging model reads

(15) \begin{equation} {\boldsymbol{y}}({\boldsymbol{x}})={\boldsymbol{F}}({\boldsymbol{x}})=f({\boldsymbol{x}})+Z({\boldsymbol{x}}) \end{equation}

The regression model $f({\boldsymbol{x}})$ can be a known or unknown constant or a multivariate polynomial, yielding various categories of the Kriging model. The correlation matrix $\boldsymbol{\Psi }$ for observation data for the Gaussian process model is defined as

(16) \begin{equation}\boldsymbol{\Psi }=\left(\begin{array}{ccc} \psi ({\boldsymbol{x}}^{(1)},{\boldsymbol{x}}^{(1)}) & \cdots & \psi ({\boldsymbol{x}}^{(1)},{\boldsymbol{x}}^{(n)})\\[5pt] \vdots & \ddots & \vdots \\[5pt] \psi ({\boldsymbol{x}}^{(n)},{\boldsymbol{x}}^{(1)}) & \cdots & \psi ({\boldsymbol{x}}^{(n)},{\boldsymbol{x}}^{(n)}) \end{array}\right)\end{equation}

where $\psi$ is the correlation function for observation data, and the most widely used correlation functions include the Gaussian, exponential, spline, linear, and spherical functions. Furthermore, the hyperparameters (including $\sigma ^{2}$ above) in the Kriging model are computed through the maximum likelihood estimation (MLE) [Reference Dietrich and Osborne29]. The predicted mean of Kriging interpolation is given by

(17) \begin{equation} \hat{{\boldsymbol{y}}}({\boldsymbol{x}})={\boldsymbol{M}}\boldsymbol{\alpha }+r({\boldsymbol{x}})\Psi ^{-1}({\boldsymbol{Y}}-\boldsymbol{\xi \alpha }) \end{equation}

where

(18) \begin{equation} {\boldsymbol{M}}=(\!\begin{array}{cccc} b_{1}(x) & b_{2}(x) & \cdots & b_{{n_{p}}}(x) \end{array}) \end{equation}

(19) \begin{equation} \boldsymbol{\alpha }=(\boldsymbol{\xi }^{T}\boldsymbol{\Psi }^{-1}\boldsymbol{\xi })^{-1}\boldsymbol{\xi }^{T}\boldsymbol{\Psi }^{-1}{\boldsymbol{Y}} \end{equation}

(20) \begin{equation} r({\boldsymbol{x}})=[\begin{array}{ccc} \psi ({\boldsymbol{x}},{\boldsymbol{x}}^{(1)}) & \cdots & \psi ({\boldsymbol{x}},{\boldsymbol{x}}^{(n)}) \end{array}] \end{equation}

and b _j is the jth polynomial regressor and α _j is the corresponding coefficient, and $0\leq j\leq n_{p}$ . $\boldsymbol{\xi }$ is the observation matrix as described above, and the estimated mean-squared error by the predictor is

(21) \begin{equation} s^{2}\!\left({\boldsymbol{x}}\right)=\sigma ^{2}\left(1-r\left({\boldsymbol{x}}\right)\boldsymbol{\Psi }^{-1}r\left({\boldsymbol{x}}\right)^{T}+\frac{1-\boldsymbol{\xi }\boldsymbol{\Psi }^{-1}r\left({\boldsymbol{x}}\right)^{T}}{\boldsymbol{\xi }^{T}\boldsymbol{\Psi }^{-1}\boldsymbol{\xi }}\right) \end{equation}

In this work, the DACE (design and analysis of computer experiments) toolbox in MATLAB [Reference Lophaven, Nielsen and Søndergaard30] is adopted for constructing the kriging model, and the linear regression model and Gaussian correlation model are used. Three Kriging models are developed for each component in the world coordinate values (X _w , Y _w , Z _w ).

2.2.5. Kriging regression model

The predictor in Eq. (21) is for Kriging interpolation, well-suited for modeling more deterministic data. In contrast, the Kriging model can also be modified for regression to model data with significant noises and uncertainties that are governed by another Gaussian process $\eta ({\boldsymbol{x}})$ with the zero mean and covariance matrix $\boldsymbol{\Sigma }$

(22) \begin{equation} \eta \sim GP(0,\boldsymbol{\Sigma }) \end{equation}

Following a similar procedure above, the corresponding Kriging regression predictor is [Reference Sobester, Forrester and Keane26]

(23) \begin{equation} \hat{{\boldsymbol{y}}}\left({\boldsymbol{x}}\right)={\boldsymbol{M}}\boldsymbol{\alpha }+r\left({\boldsymbol{x}}\right)\left(\boldsymbol{\Psi }+\frac{1}{\sigma ^{2}}\boldsymbol{\Sigma }\right)^{-1}\left({\boldsymbol{Y}}-\boldsymbol{\xi \alpha }\right) \end{equation}

The type and distribution of noises produce different covariance matrices. The simplest form of the covariance matrix describing the noise and uncertainty is

(24) \begin{equation}\boldsymbol{\Sigma }=\left(\begin{array}{ccc} \mathrm{var}(y^{(1)}) & 0 & 0\\[5pt] 0 & \ddots & 0\\[5pt] 0 & 0 & \mathrm{var}(y^{(n)}) \end{array}\right)\end{equation}

where $var({\boldsymbol{y}}^{(i)})$ is the variance of the ith data observation. Furthermore, the noise can be assumed homogeneously distributed across all the input observations, and therefore, the covariance matrix can be set as $\boldsymbol{\Sigma }=10^{\epsilon }{\boldsymbol{I}}_{n}$ , where $\epsilon$ is used to quantify the amount of noise and can be determined as a hyperparameter using the MLE above. In this study, the ooDACE Toolbox is utilized to construct the model [Reference Couckuyt, Dhaene and Demeester31]. Likewise, three kriging regression models are constructed for predicting the world coordinate X _w , Y _w , and Z _w , respectively.

4.1. 3D positioning

The model performances are first evaluated for 3D positioning and the accuracy metric is defined as [Reference Amsters, Demeester, Stevens, Lauwers and Slaets2, Reference Hong, Ou and Wang33]

(26) \begin{equation} \varepsilon =\sqrt{(X_{w,g}-\hat{X}_{w})^{2}+(Y_{w,g}-\hat{Y}_{w})^{2}+(Z_{w,g}-\hat{Z}_{w})^{2}} \end{equation}

where (X _w,g , Y _w,g , Z _w,g ) are the ground truth value of the world coordinate of the ArUco marker and are manually measured using Eq. (25) and following the procedure above. $(\hat{X}_{w},\hat{Y}_{w},\hat{Z}_{w})$ are those estimated by our models and algorithms.

Table I compares the performance of various models for all three cameras (red, white, and black, as shown in Fig. 5) in mean accuracy, lower bound, upper bound, and standard deviation. Differences in performance among the three cameras are caused by the discrepancy in camera installation and intrinsic parameters of cameras. In total, nine different models/methods are listed in Table I. Method 1 is the rigid transformation method and is denoted RT hereafter. The mean accuracy of RT for red, white, and black cameras is 1.430, 2.344, and 3.422 cm, respectively. Method 2 is the polynomial regression (PR) method, and its mean accuracy for all three cameras is around 1.5 cm. The lower bound and standard deviation of PR are both lower than 1 cm, and the upper bound is less than 4 cm. Compared to RT, PR shows excellent improvement in lower bound, upper bound, and standard deviation. In Method 3, that is, RT + PR hybrid model, RT is applied to estimate the world coordinates $(\tilde{X}_{w},\tilde{Y}_{w},\tilde{Z}_{w})$ first, and then these intermediate values are used by PR as the input to estimate the world coordinates again to further improve the accuracy. The mean accuracy, lower bound, upper bound, and standard deviation of RT + PR are the same as that of PR. Method 4, that is, ANN, is the machine learning-based regression model to estimate the world coordinates. The mean accuracy of ANN ranges from 1.5 to 2.0 cm, which is slightly worse than those of PR. The lower bound of ANN is also lower than 1 cm, and its upper bound and standard deviation are slightly higher than those of PR. In Method 5, that is, RT + ANN, the ANN model uses the estimated world coordinates from RT as inputs to estimate the world coordinates again. It seems that this hybrid model does not contribute to significant improvement in the mean accuracy. The accuracy of RT + ANN is between 1.5 and 1.8 cm, and its lower bound is also less than 1 cm. However, its upper bound and standard deviation are slightly inferior to that of ANN only. In Method 6, the Kriging Interpolation (KI) model is used to estimate the world coordinates. Compared to that of PR, the mean accuracy of KI is superior for the red camera but much worse for the white and black cameras. The upper bound and standard deviation of KI are also much higher than that of PR. Method 7 is the hybrid model combining RT + KI. Its mean accuracy is almost the same as that of KI. The differences in the lower bound, upper bound, and standard deviation between RT + KI and KI are indeed minor. Method 8 uses Kriging regression (KR) to estimate the coordinate value of the marker in the world frame. In contrast to KI, KR achieves better mean positioning accuracy for the black camera. However, the upper bound and standard deviation for the red and black cameras of KR are worse than that of KI. In Method 9, RT + KR forms a hybrid model to further improve the world coordinate estimation using the intermediate results from RT. Basically, there is no notable difference in the mean accuracy between KR and RT + KR. Their lower bound, upper bound, and standard deviation are also almost identical.

Table I. Performance of 9 different methods in 3D indoor positioning.

Through the comparison of these data-driven modeling methods, several interesting observations can be made. The combination of RT and other methods, such as RT + PR, RT + ANN, RT + KI, and RT + KR, makes a negligible contribution to mean accuracy improvement when compared to the single PR, ANN, KI, and KR. The difference in the mean accuracy among these nine methods is minor for the red camera but is noticeable for white and black cameras. PR and RT + PR are clearly top performers and achieve a mean accuracy of around 1.5 cm for all cameras and greatly outperform the other methods. And the lower bound, upper bound, and standard deviation for PR and RT + PR are the same. Compared to existing research works using the ArUco marker, the present system achieves a value of around 1.5 cm, a significant improvement in 3D indoor positioning accuracy.

All these nine modeling methods in Table I perform well in positioning when the number of training data is 40 for each camera. However, the manual measurement and data collection are tedious and time-consuming and should be minimized subject to the positioning accuracy requirement. Therefore, the tradeoff between positioning accuracy and the training data size is also studied. As shown in Fig. 10, an analysis is conducted to investigate the influence of the number of training data points (10, 20, 30, and 40) on positioning accuracy. Figure 10(a), (b), and (c) show the mean accuracy for the red, white, and black cameras, respectively. Again, the training data points are selected randomly from the data pool in Fig. 9, and the locations of testing points remain the same. It is interesting to observe that the mean accuracy of these methods, in general, improves as more training data points are incorporated. All of them yield a mean accuracy of less than 10 cm when the number of training data points is 20 or more. In addition, even with the same method, the accuracy of different cameras can be different because the locations of selected training and testing data and the camera installation all affect model performance. The mean accuracy of ANN, RT + ANN, KR, and RT + KR deteriorates appreciably when training data are limited, while that of RT, PR, RT + PR, KI, and RT + KI maintain relatively constant with different amounts of training data. Especially, RT performs very consistently and even exceeds PR when only 10 data points are used, and its mean accuracy for all three cameras mostly keeps below 5 cm, which may be ascribed to physics/kinetics contained in its mathematical form in Eq. (6) that alleviates the demand for data. Our observation implies that when training data is scarce, and it is difficult to collect more, RT and KI may be more appropriate for modeling the mapping relationship F .

Figure 10. The effect of the number of training data points on 3D positioning accuracy: (a) red, (b) white, and (c) black cameras.

4.2. 2D positioning

When the system is used for a mobile robot, and the height of the attached ArUco marker does not change appreciably, the estimation of the Z _w coordinate value is not necessary. Thus, the 2D positioning accuracy is more important, and its evaluation metric will be rewritten by removing the term associated with Z _w in Eq. (26)

(27) \begin{equation} \varepsilon =\sqrt{(X_{w,g}-\hat{X}_{w})^{2}+(Y_{w,g}-\hat{Y}_{w})^{2}} \end{equation}

In Table II, the mean accuracy, lower bound, upper bound, and standard deviation for the nine methods in the 2D domain are listed. The mean accuracy of RT varies between 0.9 and 1.8 cm. The lower bound of RT is relatively small and only 0.4 cm. However, the upper bound is large and reaches 7.126 cm. Besides, its standard deviation is also larger than most methods, as shown in Table II. Among all methods, PR and RT + PR again perform the best and achieve the 2D mean accuracy of about 1 cm, and their lower bound, upper bound, and standard deviation are also smaller than those of the other methods. The mean accuracy of ANN, RT + ANN, KI, and RT + KI ranges from 1.1 to 1.8 cm, which is somewhat worse than that of PR and RT + PR. Correspondingly, their lower bound, upper bound, and standard deviation are slightly worse relative to PR and RT + PR. The mean accuracy of KR and RT + KR is the worst compared to other methods, and the discrepancy between the estimation and ground truth can reach 2.307 cm. Although their lower bound is lower, the upper bound and the standard deviation are high and, respectively, over 8 and 2.5 cm. Again, PR and RT + PR exhibit excellent positioning performance, and the combination of RT and PR does not offer apparent advantages over the PR only.

Table II. Performance of 9 different methods in 2D indoor positioning.

The effect of the number of training data on the mean accuracy of 2D positioning for red, white, and black cameras is shown in Fig. 11(a), (b), and (c), and the general trend and dependence on the training data size is similar to that of 3D positioning. The methods of ANN, RT + ANN, RK, and RT + RK perform poorly when the number of training data points is only 10. All these methods perform very well with a mean accuracy smaller than 10 cm when 20 or more data points are used for model training. RT, RT + PR, PR, KI, and RT + KI are able to keep relatively consistent positioning performance even with only ten training data points, while RT is particularly appealing. Similarly, when only a limited amount of data is available for model training, RT, KI, or RT + KI are preferred.

Figure 11. The effect of the number of training data points on 2D positioning accuracy: (a) red, (b) white, and (c) black cameras.

4.3. Verification with mobile robot experiment

Next, a test is carried out to verify that the indoor positioning system can be used to track moving robots in real time for research uses. The indoor positioning system proposed in this study is implemented on a laptop with an Intel(R) Core (TM) i5-7200U CPU @ 2.50 GHz. The mobile robot moves at a maximum speed of 0.26 m/s. The system processes three frames simultaneously from the three cameras within approximately 0.05 s, which allows real-time positioning and control of the robot. The ArUco marker is placed on the top of a mobile robot, and the centers of both are aligned. As shown in Fig. 12, “+” labels are placed at four different locations, which are, respectively, P ₁(0.5, 0.5), P ₂(0.5, 2.5), P ₃(2.5, 2.5), and P ₄(4.5, 4.5), all with a unit of meter. They are set as the waypoints for the mobile robot, and the starting position of the robot is around P ₁(0.5, 0.5), as shown in Fig. 12(a). In addition, the estimated location of the marker’s center in the world coordinate frame, that is, $\hat{X}_{w0}$ and $\hat{Y}_{w0}$ is shown at the top-left corner of each subfigure in Fig. 12. Then the onboard PID controller uses the position information of $\hat{X}_{w0}$ and $\hat{Y}_{w0}$ to command the robot to go through the four waypoints above sequentially. In Fig. 12(b), it is clearly seen that the mobile robot moves from P ₁(0.5, 0.5) to P ₂(0.5, 2.5). In Fig. 12(c) and (d), the robot reaches P ₂(0.5, 2.5) and P ₃(2.5, 2.5) in tandem. Subsequently, the robot is on its way toward P ₄(4.5, 4.5), as shown in Fig. 12(e). Finally, the robot arrives at destination P ₄(4.5, 4.5) in Fig. 12(f). It should be noted that because of the partially overlapped view between two adjacent cameras, P ₂ is present in both the white (middle) and black (right) cameras, and sometimes the robot appears simultaneously in two camera images (e.g., Fig. 12(b), (c), and (e)). This experiment verifies that the present indoor positioning system can track the mobile ground robot and provides accurate position information in real time for robot control.

Figure 12. Experiment to verify the feasibility of using our indoor positioning system for robot control and research uses.