Accurate and stable alignment of virtual and real spaces using consumer-grade trackers

Abstract

The world space defined by a traditional virtual reality application is generally regarded as separate from the real-world space in which users actually exist. In this paper, we present a method that enables such detached spaces to connect in an integrated space. In particular, we show that by using a specially manufactured calibration board and three consumer-grade position tracking devices, a reference coordinate system can easily be set up in the physical space, whose geometric relationship with the virtual reality space is estimated with high numerical accuracy and stability. Then, we demonstrate that, combined with traditional computer vision techniques for marker tracking, the presented technique allows colocated users from virtual, augmented, and mixed realities to cooperate with each other while making effective use of technologies from other realities.

Appendices

Appendix 1: The numerical method to solve the nonlinear least-squares problem in Sect. 3.2.2

Consider a triangle (not necessarily equilateral) on the plane formed by three 2D points $\tilde{\mathbf{p}}_i = ({\tilde{x}}_i, {\tilde{y}}_i)^{\top } (i = 0, 1, 2)$, whose center of mass is assumed to exist, without loss of generality, at the origin of the 2D coordinate system (see Fig. 4 again). Now, the problem is to find an equilateral triangle of side $s = 200\sqrt{3}$ mm that fits best, in the least-squares sense, the given triangle. This can be done by computing a rigid body transformation that transforms the triangle, defined by three vertices $\mathbf{q}_i = (x_i, y_i)^{\top }$ with $\mathbf{q}_0 = (\frac{1}{2}s, -\frac{\sqrt{3}}{6}s)^{\top }$, $\mathbf{q}_1 = (0, \frac{\sqrt{3}}{3}s)^{\top }$, and $\mathbf{q}_2 = (-\frac{1}{2}s, -\frac{\sqrt{3}}{6}s)^{\top }$, so that the sum of squared distances between the corresponding points is minimized.

Let $M(\alpha , \beta ;\theta ) = T(\alpha , \beta )\cdot R(\theta )$ be a 2D rigid body transformation represented by a rotation followed by a translation. Then, for the transformed positions $\bar{\mathbf{q}}_{i}$ of $\mathbf{q}_{i}$ with $\bar{\mathbf{q}}_{i} = (x_{i} \cos \theta - y_{i} \sin \theta + \alpha , \, x_{i} \sin \theta + y_{i} \cos \theta + \beta )^{\top }$, the squared sum of the distances between $\bar{\mathbf{q}}_{i}$ and $\bar{\mathbf{p}}$ is expressed as

$$\begin{aligned} F(\alpha , \beta ;\theta )& {}= \sum _{i=0}^{2} \big \{ ( x_{i} \cos \theta - y_{i} \sin \theta + \alpha - {\tilde{x}}_{i})^{2} \nonumber \\&\quad + (x_{i} \sin \theta + y_{i} \cos \theta + \beta - {\tilde{y}}_{i})^{2} \, \big \}. \end{aligned}$$

Then, by noting that $\sum _0^2 x_i = \sum _0^2 y_i = 0$, we are led to

$$\begin{aligned} \frac{\partial F}{\partial \alpha } = \sum _{i=0}^{2} {2 (x_{i} \cos \theta - y_{i} \sin \theta + \alpha - {\tilde{x}}_{i}) } = 2\left (3\alpha - \sum _{i=0}^{2} {\tilde{x}}_{i}\right). \end{aligned}$$

From $\nabla F(\alpha , \beta ;\theta ) = 0$ and the similarity between the two variables $\alpha$ and $\beta$, we find the summed error is minimized at $({\alpha }^*, {\beta }^*)$ where ${\alpha }^* = \sum _{i=0}^{2} {\tilde{x}}_{i}/{3}$ and ${\beta }^* = \sum _{i=0}^{2} {\tilde{y}}_{i}/{3}$. For the other variable $\theta$,

$$\begin{aligned} \frac{\partial F}{\partial \theta } = 2 \left \{ \,\sin \theta \sum _{i=0}^{2} ({\tilde{x}}_{i}x_i + {\tilde{y}}_{i}y_{i}) + \cos \theta \sum _{i=0}^{2} ({\tilde{x}}_{i}y_i - {\tilde{y}}_{i}x_{i}) \,\right \}. \end{aligned}$$

To obtain ${\theta }^*$ such that $\frac{\partial F}{\partial \theta } = 0$, the Newton–Raphson method (Solomon 2015) is applied to the nonlinear function $f(\theta ) = S\cdot \sin \theta + C\cdot \cos \theta$ with $S = \sum _{i=0}^{2} ({\tilde{x}}_{i}x_i + {\tilde{y}}_{i}y_{i})$ and $C = \sum _{i=0}^{2} ({\tilde{x}}_{i}y_i - {\tilde{y}}_{i}x_{i})$. Note that, since the rotational displacement is usually very small, the Newton–Raphson iteration converges quickly to the solution ${\theta }^*$ for initial value $\theta _0 = 0$:

$$\begin{aligned} \theta _{k+1} = \theta _{k}-\frac{f(\theta _{k})}{f'(\theta _{k})} = \theta _{k}-\frac{S \sin \theta _{k} + C \cos \theta _{k}}{S \cos \theta _{k} - C \sin \theta _{k}}. \end{aligned}$$

Appendix 2: The Gauss–Newton method to solve the nonlinear least-squares problem in Sect. 3.2.3

Consider an enumeration of the entire position samples $\mathbf{p}_{ij}$ ($i = 0, 1, 2$ and $j = 0, 1, \ldots , n_{{\rm total}}-1$), for which i(k) and j(k) denote the tracker number i and the tracker’s sample number j of the kth sample, respectively. If we define $\mathbf{f}({\varvec{\xi }})$ to be a function from ${\mathbb {R}}^6$ to ${\mathbb {R}}^{9\cdot n_{{\rm total}}}$ whose value is a vector constructed by stacking the vectors $\tau (\mathbf{q}_{ i(k)},T({{\varvec{\xi }}})) - \mathbf{p}_{i(k),j(k)} \in {\mathbb {R}}^3$ in the enumeration order, the nonlinear minimization problem can be restated as

$$\begin{aligned} {{\varvec{\xi }}}_{{\rm opt}} = \arg \min _{{\varvec{\xi }}}\, \mathbf{f}({\varvec{\xi }})^{\top }\mathbf{f}({\varvec{\xi }}). \end{aligned}$$

Then, given an initial value ${\varvec{\xi }}^{(0)}$, the current estimate for the transformation parameter is improved iteratively until convergence as ${\varvec{\xi }}^{(k+1)} = {\varvec{\xi }}^{(k)} + \Delta \,{\varvec{\xi }}$ by solving the $6\times 6$ linear system for $\Delta {\varvec{\xi }}$

$$\begin{aligned} \mathbf{J}_\mathbf{f}({\varvec{\xi }}^{(k)})^{\top }{} \mathbf{J}_\mathbf{f}({\varvec{\xi }}^{(k)})\Delta \,{\varvec{\xi }}= -\mathbf{J}_\mathbf{f}({\varvec{\xi }}^{(k)})^{\top }{} \mathbf{f}({\varvec{\xi }}^{(k)}), \end{aligned}$$

where $\mathbf{J}_\mathbf{f}({\varvec{\xi }}^{(k)}) \in {\mathbb {R}}^{9\cdot n_{{\rm total}}\times 6}$ is the Jacobian matrix of $\mathbf{f}$ at ${\varvec{\xi }}^{(k)}$, which can be constructed by stacking the following $3\times 6$ matrix for the kth position sample (Claraco 2018):

$$\begin{aligned} \left[ \begin{array}{c|c} I_{3\times 3}&-[\tau (\mathbf{p}_{i(k),j(k)},T({\varvec{\xi }}^{(k)}) )]_{\times } \end{array} \right] . \end{aligned}$$

(Here, $I_{3\times 3}$ and $[\cdot ]_{\times }$ represent the $3\times 3$ identity matrix and the skew-symmetric matrix operator, respectively.) Note that when the parameter ${\varvec{\xi }}$ is updated in each iteration, the rigid body transformation is also updated correspondingly as $T({\varvec{\xi }}^{(k+1)}) = T(\Delta {\varvec{\xi }})T({\varvec{\xi }}^{(k)})$.