3.1.1. The human visual system and gaze tracking

Humans perceive visual information through sensory receptors in the eyes. The process begins when light passes through the cornea, enters the pupil, and then gets focused by the lens onto the retina. This is then processed in the brain where an image is formed. The retina has two kinds of photoreceptors: cones and rods. Cones are capable of color vision and are responsible for high spatial acuity. Rods are responsible for vision at low light levels. As shown in Figure 4-A, the cone density is highest in the central region of the retina and reduces monotonically to a fairly even density in the peripheral retina region. Retinal eccentricity (or simply, eccentricity) implies the angle at which the light from the image gets focused on the retina. This distribution of the photoreceptors gives rise to the concept of Foveation and helps define the idea of visual acuity.

The density of photoreceptors (cones and rods) declines monotonically and in a continuous manner from the center of the retina to its periphery. Nevertheless, approximating the retina as being formed of discrete concentric regions, where the density of the photoreceptors corresponds to eccentricity angles, helps simplify the concept – this is the concept of Foveation. A summary of photoreceptors distribution is presented in Table II, followed by further details below:

• • The region from the center of the retina up to 5 $^{\circ }$ of eccentricity is the Fovea region. The Fovea is only about 1% of the retina but has the highest density of cone photoreceptor cells, and the brain’s visual cortex dedicates about 50% of its area to information coming from the Fovea [Reference Sherman, Craig, Sherman and Craig45]. Therefore, the Fovea has the highest sensitivity to fine details.

• • The region that surrounds the Fovea is commonly known as the Parafovea, which goes up to 8 $^{\circ }$ of the visual field [Reference Hendrickson, Penfold and Provis16]. The parafoveal region provides visual information as to where the eyes should move next (saccade) and supports the Fovea to process the region of interest in detail. Previous research has investigated if meaningful linguistic information can be obtained from parafoveal visual input while reading [Reference Hyönä, Liversedge, Gilchrist and Everling46].

• • The next region that surrounds the Parafovea is called Perifovea, which extends approximately up to 18 $^{\circ }$ of eccentricity. In this region, the density of rods is higher than that of cones, about 2:1. Consequently, unlike the Fovea and Parafovea, only rough changes in shapes are perceived in this region [Reference Ishiguro and Rekimoto47].

• • The region beyond 18 $^{\circ }$ , and up to about 30 $^{\circ }$ of the visual field, is known as the Near-Peripheral Region. It has the distribution of 2–3 rods between cones [Reference Quinn, Csincsik, Flynn, Curcio, Kiss, Sadda, Hogg, Peto and Lengyel44]. This region is responsible for the segmentation of visual scenes into texture-defined boundaries (“texture segregation”) and the extraction of contours for pre-processing in pattern and object recognition [Reference Strasburger, Rentschler and Jüttner48].

• • The region between 30 $^{\circ }$ and 60 $^{\circ }$ of eccentricity is called the Mid-Peripheral Region [Reference Simpson49]. Although acuity and color perception degrade rapidly in this region, researchers have shown that color perception is still possible even at large eccentricities, up to $\sim$ 60 $^{\circ }$ [Reference Gordon and Abramov50].

• • The region at the edge of the visual field (from 60 $^\circ$ up to nearly 180 $^\circ$ horizontal diameter) is called the Far Peripheral Region. This region has widely separated ganglion cells, and visual functions such as stimulus detection, flicker sensitivity, and motion detection are still possible here [Reference Strasburger, Rentschler and Jüttner48].

Table II. Human retinal regions and their sizes in diameter and angles (derived from [Reference Quinn, Csincsik, Flynn, Curcio, Kiss, Sadda, Hogg, Peto and Lengyel44]).

Visual acuity can be quantitatively represented in terms of minimum angle of resolution ( $\mathrm{MAR}_{0}$ , measured in arcminutes) [Reference Guenter, Finch, Drucker, Tan and Snyder17, Reference Strasburger, Rentschler and Jüttner48, Reference Weymouth51]. $\mathrm{MAR}_{0}$ can be understood as the smallest angle at which two objects in the visual scene are perceived as separate [Reference Weymouth51]. $\mathrm{MAR}_{0}$ accounts for the number of neurons allocated to process the information from the visual field, as a function of the eccentricity. This relation between $\mathrm{MAR}_{0}$ and eccentricity can be approximated as a linear model, which has been shown to closely match the anatomical features of the eye [Reference Guenter, Finch, Drucker, Tan and Snyder17, Reference Bruder, Schulz, Bauer, Frey, Weiskopf and Ertl37, Reference Strasburger, Rentschler and Jüttner48, Reference Weymouth51].

(1) \begin{equation} \mathrm{MAR} = mE + \textit{MAR}_{0} \end{equation}

Here, ${\mathrm{MAR}_0}$ is the intercept, which signifies the smallest resolvable eccentricity angle for humans, and $m$ is the slope of the linear model. ${\mathrm{MAR}_0}$ for a healthy human varies between 1 and 2 arcminutes, that is, $1/60^\circ$ to $1/30^\circ$ (1 $^\circ$ = 60 arcminutes). Authors in ref. [Reference Guenter, Finch, Drucker, Tan and Snyder17] experimentally determined the values of $m$ based on observed image quality, ranging between 0.022 and 0.034. Figure 4-B captures this linear relationship of Eq. (1), showing how visual acuity degrades as a function of eccentricity, represented with a piece-wise constant approximation, that is, each retinal region has a distinct constant $\mathrm{MAR}$ value [Reference Guenter, Finch, Drucker, Tan and Snyder17].

Figure 5. A schematic showing the coordinate system on the remote and operator sites.

The information from the operator’s site is subsequently transmitted to foveate the 3D point cloud. Foveating a 3D point cloud implies introducing concentric regions in it that correspond to the retinal fovea regions. The regions are centered on the human eye-gaze direction, each of them having a specific radius and its associated visual acuity, that is, rendering quality.

3.1.2. Visualization and coordinate transformations

To visualize and explore the incoming point cloud, as well as the remote robotic platforms, a VR-based interface is designed using the Unreal graphics engine on Windows 10. This creates the immersive remote teleoperation environment for the operator. As noted earlier, there is an interdependence between the operator site and the remote site, in that, the eye-gaze data from the operator site is required for the foveation model at the remote site. The foveated point cloud is then streamed back to the operator site to be rendered for visualization. Furthermore, the operator site and the remote site are independent environments with their respective reference frames. It is therefore necessary to implement appropriate transformations among all the entities to ensure correct data exchange and conversion. As shown in Figure 5, the reference frames are as follows: UE world coordinate frame $\textbf{U}$ , HMD coordinate frame $\textbf{E}$ , and the gaze direction vector $^E\vec{D}$ $\in R^3$ on $\textbf{E}$ .

To calculate the correct gaze pose in $\textbf{U}$ , transforming $^E\vec{D}$ from $\textbf{E} \to \textbf{U}$ is required, through the head pose $^U\textbf{H}$ on U, as follows:

(2) \begin{equation} ^{U}\vec{D} = ^U\textbf{H}\cdot ^E\vec{D} \end{equation}

This gaze direction $^{U}\vec{D}$ , along with the head pose $\textbf{H}$ are communicated to the remote site for additional processing. Further, the received point cloud from the remote site is visualized in Unreal and needs to be positioned based on the pose of the camera positioned at the remote site. At the remote site, the coordinate system of the camera pose, $^O\mathbf{P}$ is in OpenGL, $\textbf{O}$ . The pose has to be transformed, using a change-of-basis matrix, into the UE coordinate system. UE uses a left-handed, $z$ -up coordinate system, while the camera coordinates of OpenGL use a right-handed coordinate system, $y$ -up. Eq. (3) provides the coordinate transformation formula, where $\mathbf{B}$ is the change-of-basis transformation matrix.

(3) \begin{equation} ^U\mathbf{P} = \mathbf{B}\cdot ^O\mathbf{P}\cdot \mathbf{B}^{-1} \end{equation}

At the operator site, rendering the received real-time dynamic point-cloud data from the remote site requires a high-speed large data transfer, as well as efficient and high-quality visualization. To meet these requirements, the following modules were developed, as seen in Figure 3:

A real-time point-cloud decoder: that decompresses the data received at the operator site. The decoding module includes the state-of-the-art point-cloud codec algorithm from ref. [Reference Mekuria, Blom and Cesar26] and implemented the Boost ASIO over a TCP socket for data transfer. As the point cloud codec needs to compress and stream point clouds using gaze information from master station, TCP was chosen to ensure data integrity and packet reception.

Conversion system: Each decoded point-cloud region $\boldsymbol{\mathcal{P}}_n (\forall{n} \in{1,\cdots \,,N})$ has to be converted into a texture for visualization, where the reference frame of the received data has to be transformed into that of the user site, that is, the unreal graphics engine coordinate system.

A rendering system: The data should be transferred to the GPU and made accessible to the graphics engine shader for real-time rendering. We’ve implemented a splatting-based technique to render point clouds using UE’s Niagara particle system. After decoding the point cloud data, positions, and colors, we transferred this data to the GPU via the Niagara module. Once on the GPU, each particle’s UV coordinates are calculated based on its unique identifier, texture size, and virtual camera positions. This process ensures precise mapping of the point cloud data onto the rendered surfaces. By integrating Niagara’s rendering capabilities, we achieve faster rendering performance, which is around eight milliseconds.

3.1.3. Robot control

The motion controllers of the HTC Vive Pro Eye allow remote teleoperators to send control commands in real-time using a wireless motion controller interface. These HTC Motion Controllers (HMC) are tracked in space, allowing the teleoperator to freely explore the VR environment and interact naturally with the remote robots. The HMC has two buttons: To engage/disengage motion commands between the HMC and remote robot, overcoming range limitations, and to open/close the remote robot’s end-effector for precise object manipulation and control. Since the operator directly commands the motion of the remote robot using a HMC controller, our design follows a human-in-the-loop design paradigm without assuming any autonomy or semi-autonomy of the robots [Reference Naceri, Mazzanti, Bimbo, Tefera, Prattichizzo, Caldwell, Mattos and Deshpande4]. The remote environment point cloud is sampled, streamed, and rendered for the human operator in real-time. Based on the remote scene, the operator plans the next step and the action commands are sent to the remote robot. The remote robots are considered passive and do not act autonomously. Velocity controller was used due to its ability to provide smooth and continuous control over the movement of a remote robotic system. It calculates joint angles, velocities, and accelerations, which are then transmitted from the HMC interface to the remote robots. The joint states of the remote robots are relayed back to the interface to update the poses of the virtual robot’s models. This integration allows the operator to see and command the motion of the virtual robot, with the motion of the real robot being mapped 1:1 to the virtual robot motion. Due to the non-homothetic nature of the kinematics between the remote robot and the operator controllers, the velocity-control commands the 7-DOF pose of the robot. The inverse Jacobian method is employed to calculate the velocity of the HMC, which is then mapped to the robot. To command the robot pose, the position error $\mathbf{e} \in SE(3)$ is determined by comparing the current pose with the desired pose. The damped least-squares solution, as shown in Eq. (4), is iteratively used to find the change in joint angles $\Delta \mathbf{q}$ that minimizes the error $\mathbf{e}$ .

(4) \begin{equation} \Delta \mathbf{q} = (\mathbf{J}^T\mathbf{J} + \lambda \mathbf{I})^{-1}\mathbf{J}^T \mathbf{e}, \end{equation}

where $\mathbf{J}$ is the manipulator Jacobian, $\lambda \in R$ is a non-zero damping constant, and $\mathbf{I}$ is the identity matrix. Finally, a proportional controller $\dot{\mathbf{q}} = K_p \cdot \Delta \mathbf{q}$ sets the joint velocities to achieve the desired pose. The values for $\lambda$ and $K_p$ were determined empirically as 0.001 and 0.6, respectively.

3.2.1. Coordinate transformations

As shown in Figure 3, a parallel module will receive information about head pose $^U\mathbf{H}$ and the gaze direction $^U\vec{D}$ from the Operator site. The coordinate system of the head pose $^U\mathbf{H}$ and the gaze direction $^U\vec{D}$ has to be transformed from the Unreal frame, $\textbf{U}$ , to the remote site OpenGL coordinate system, $\textbf{O}$ , using a similar change-of-basis matrix, Eq. (5). Here: $\mathbf{Q}$ is the change-of-basis transformation matrix from UE to OpenGL. Figure 5 left shows the reference frames of the remote site, and they are described as follows: OpenGL world coordinate frame $\textbf{O}$ and the camera frame $\textbf{C}$ .

(5) \begin{equation} \begin{split} ^O\mathbf{H} = \mathbf{Q}\cdot ^U\mathbf{H}\cdot \mathbf{Q}^{-1} \\ ^O\vec{D} = \mathbf{Q}\cdot ^U\vec{D}\cdot \mathbf{Q}^{-1} \end{split} \end{equation}

Gaze direction vector $^O\vec{D}$ and the head pose $^O\mathbf{H}$ have to be transformed into the camera coordinate frame, in order to perform 3D reconstruction, map partitioning, and sampling. For every frame, the color image $\mathbf{C}$ and depth map $\mathbf{D}$ are registered into the map model $\boldsymbol{\mathcal{M}}$ by estimating the global pose of the camera $\textbf{P}$ . Section 3.2.2 will provide a concise description of this.

(6) \begin{equation} \begin{split} ^{C}\mathbf{H} = \textbf{P}^{-1}\cdot ^O\mathbf{H}\\ ^C\vec{D} = \textbf{P}^{-1}\cdot ^O\vec{D} \end{split} \end{equation}

The head pose in the camera frame is used as a point of gaze origin H (hx,hy,hz) and using the gaze direction vector, $^C\vec{D}$ , a ray, that is, the gaze vector $\textbf{L}$ is projected into the 3D map.

3.2.2. 3D data acquisition, mapping, partitioning, and sampling

The acquisition and point generation module acquires RGB-D images from the RGB-D cameras, for example, Intel RealSense and ZED stereo camera. The generated points (maps) pipeline leverages the state-of-the-art real-time dense visual SLAM system, ElasticFusion [Reference Whelan, Leutenegger, Salas-Moreno, Glocker and Davison24] for initial camera tracking. The map $\boldsymbol{\mathcal{M}}$ is represented using an unordered list of surfels, where each surfel $\boldsymbol{\mathcal{M}}^s$ has a position $\mathbf{p} \in \mathbb{R}^3$ , a normal $\mathbf{n} \in \mathbb{R}^3$ , a color $\textbf{c}$ $\in$ $\mathbb{R}^3$ , a weight $w$ $\in$ $\mathbb{R}$ , a radius $r$ $\in$ $\mathbb{R}$ , an initialization timestamp $t_0$ , and a current timestamp $t$ . The camera intrinsic matrix K is defined by: (i) the focal lengths $f_x$ and $f_y$ in the direction of the camera’s $x-$ and $y-$ axes, (ii) a principal point in the image $(c_x,c_y)$ , and (iii) the radial and tangential distortion coefficients $k_1,k_2$ and $p_1,p_2$ respectively. The domain of the image space in the incoming RGB-D frame is defined as $\Omega$ $\subset{\mathbb{N}}^2$ , with the color image $\mathbf{C}$ having pixel color c : $\Omega \to \mathbb{N}^3$ , and the depth map $\mathbf{D}$ having pixel depth $d$ : $\Omega \to \mathbb{R}$ .

A. 3D point partitioning: For brevity, the symbol $\boldsymbol{\mathcal{M}}$ is used for real-time point-cloud. The density of $\boldsymbol{\mathcal{M}}$ , especially at high resolutions, implies increased computational complexity and more graphical and time resources for streaming it in immersive remote teleoperation. The foveation model can be utilized here to reduce the data. By projecting the retinal fovea regions into it, $\boldsymbol{\mathcal{M}}$ is partitioned into regions. It is then resampled to approximate the monotonically decreasing visual acuity in the foveation model, termed foveated sampling.

To partition $\boldsymbol{\mathcal{M}}$ into regions, the discussion in Section 3.1.1 is taken forward. The center of the eye gaze is used as a point of origin H $(hx,hy,hz)$ . To partition $\boldsymbol{\mathcal{M}}$ into $\boldsymbol{\mathcal{M}}_n$ regions $\forall{n} \in \{0\cdots \,N \}$ , for each of the $N$ retinal regions, a ray is cast from H $(hx,hy,hz)$ . This ray, i.e., the gaze vector $\textbf{L}\in{\mathbb R}^3$ is extended up to last point of intersection G $(gx,gy,gz)$ with the surfel map. The foveated regions are now structured around L and (Figure 1) shows the foveated regions sampled, streamed and rendered. With 3D data, the concentric regions are conical volumes, with their apex at H $(hx,hy,hz)$ , with increasing radii away from H. Algorithm1, which is implemented in Compute Unified Device Architecture (CUDA) in the GPU for faster processing, details how the radii are calculated based on $d^{vi}$ for each surfel. To assign each surfel in $\boldsymbol{\mathcal{M}}$ to a particular region $\boldsymbol{\mathcal{M}}_n$ , the shortest distance, i.e., the perpendicular distance between the surfel and L is used. As shown in Figure 1, the shortest distance from the surfel P $(px,py,pz)$ to the ray $\textbf{L}$ is the perpendicular $\textbf{PB} \bot \textbf{L}$ , where B $(bx,by,bz)$ is a point on $\textbf{L}$ . $\lVert PB\rVert$ can be obtained using the the projection of $\vec{HP}$ on $\textbf{L}$ , that is, the cross product of $\vec{HP}$ and $\vec{HG}$ , normalized to the length of $\vec{HG}$ , refer Eq. (7). Algorithm1 assigns surfel $\textbf{P}$ to the region $\boldsymbol{\mathcal{M}}_n$ .

(7) \begin{equation} \lVert PB\rVert _{\textbf{P}} = \frac{ \lVert{\vec{HP} \times \vec{HG}}\rVert }{\lVert{\vec{HG}}\rVert } \end{equation}

Algorithm 1. 3d Point Partitioning Algorithm

B. Foveated Point-Cloud (PCL) Sampling: The partitioned global map $\boldsymbol{\mathcal{M}}$ , with the region-assigned surfels, is converted into a PCL point-cloud data structure, $\boldsymbol{\mathcal{P}}_n$ for each $\boldsymbol{\mathcal{M}}_n$ region $\forall{n} \in \{ 0\cdots \,N \}$ . This conversion is sped up using a CUDA implementation in the GPU. To implement the foveated sampling, the $\mathbb{R}^3$ space of each $\boldsymbol{\mathcal{P}}_n$ region needs to be further partitioned into an axis-aligned regular grid of cubes. This process of re-partitioning the regions is called voxelization and the discrete grid elements are called voxels. After voxelization, the down-sampling of the PCL follows the foveation model – the voxels in the fovea region of the PCL are the densest, and this density progressively reduces toward the peripheral regions.

This voxelization and down-sampling is a three-step process: (1) calculating the volume of the voxel grid in each region, (2) calculating the voxel size, that is, dimension, $\mathfrak{v}_n$ , for the voxelization in each region, and (3) down-sampling the point cloud inside each voxel for the region by approximating it with the 3D centroid point of the point cloud.

Calculating the volume of the voxel grid for each region is done by simply calculating the point cloud distribution for that region, $[(x_{n,min},x_{n,max}),(y_{n,min},y_{n,max}),(z_{n,min},z_{n,max})]$ . Calculating the voxel size, $\mathfrak{v}$ , is a more involved process. Here, the visual acuity discussion from Section 3.1.1 is utilized. Consider the voxelization of the central fovea region, $\boldsymbol{\mathcal{P}}_0$ . As noted earlier, the smallest angle a healthy human with a normal visual acuity of 20/20 can discern is 1 arcminute, that is, 0.016667 $^\circ$ . Following Eq. (1) therefore, $\mathrm{MAR}_0$ , which is the smallest resolvable angle, is $0.016667^\circ$ . The smallest resolvable object length on a virtual image can be calculated as:

(8) \begin{equation} \mathfrak{l} = d^{vi} * \tan \left (\mathrm{MAR}_0\right ) \end{equation}

Equation (8) itself could provide the optimum voxel size, $\mathfrak{v}$ . The important consideration here is the value of $d^{vi}$ . In Algorithm 1, a $d^{vi}$ value for each point is calculated. In contrast, here in order to down-sample the region based on the voxelization, we calculate one $d^{vi}$ value for the entire $\boldsymbol{\mathcal{P}}_0$ region, approximated as the distance from the eye (point of gaze origin) to the 3D centroid of the point-cloud in the region, Eq. (9).

(9) \begin{equation} \mathfrak{pc}_{0} = \frac{1}{N_{\boldsymbol{\mathcal{P}}_0}}\left ({\sum _{i=1}^{N_{\boldsymbol{\mathcal{P}}_0}}{x_i}},{\sum _{i=1}^{N_{\boldsymbol{\mathcal{P}}_0}}{y_i}},{\sum _{i=1}^{N_{\boldsymbol{\mathcal{P}}_0}}{z_i}} \right ) \end{equation}

(10) \begin{equation} d_0^{vi} = \mathbf{d}(\textbf{H}, \mathfrak{pc}_{0}) \end{equation}

where $N_{\boldsymbol{\mathcal{P}}_0}$ is the number of PCL points in $\boldsymbol{\mathcal{P}}_0$ , and $\textbf{H}$ is the eye-gaze origin. Then, Eq. (8) is re-written as Eq. (11) to give the voxel size $\mathfrak{v}_0$ for the region.

(11) \begin{equation} \mathfrak{v}_0 = d_0^{vi} * \tan \left (\mathrm{MAR}_0\right ) \end{equation}

Figure 6. Foveated point-cloud sampling example showing the different voxel grid sizes for the different regions.

Once the voxelization of region $\boldsymbol{\mathcal{P}}_0$ is finalized, for the subsequent concentric regions from $\boldsymbol{\mathcal{P}}_1$ to $\boldsymbol{\mathcal{P}}_n$ , the voxel sizes are correlated with the linear $\mathrm{MAR}$ relationship in Figure 4. Eq. (12) shows that as the eccentricity angle of the regions increases, so do the voxel sizes.

(12) \begin{equation} \begin{split} \mathrm{MAR}_{n} = m\cdot E_{n} + \mathrm{MAR}_0 \\[5pt] \mathfrak{v}_{n} = \frac{\mathrm{MAR}_n}{MAR_{n-1}} * \mathfrak{v}_{n-1} \end{split} \end{equation}

The increasing voxel size away from the fovea region implies more and more points of the point cloud of the corresponding regions are now accommodated within each voxel of that region. Therefore, when the down-sampling step is applied, the approximation of the point cloud within a voxel is done over progressively dense voxels. For the down-sampling part, the region $\boldsymbol{\mathcal{P}}_0$ being the fovea region is left untouched so its density is the same as the incoming global map density. The down-sampling in the subsequent regions is done by approximating the point cloud within each voxel with its 3D centroid, using Eq. (13).

(13) \begin{equation} \mathfrak{pc}_{n}^v\left (x,y,z\right ) = \frac{1}{N_{\boldsymbol{\mathcal{P}}_{n}}^v}\left ({\sum _{i=1}^{N_{\boldsymbol{\mathcal{P}}_{n}}^v}{x_i}},{\sum _{i=1}^{N_{\boldsymbol{\mathcal{P}}_{n}}^v}{y_i}},{\sum _{i=1}^{N_{\boldsymbol{\mathcal{P}}_{n}}^v}{z_i}} \right ) \end{equation}

Here, $N_{\boldsymbol{\mathcal{P}}_{n}}^v$ is the number of points in voxel $v$ of the region $\boldsymbol{\mathcal{P}}_{n}$ $\left (\forall{n} \in \{{1\cdots \,N}\}\right )$ . Figure 6 shows the sample voxel grids for the different regions. The foveated sampling and compression is the most computationally expensive system component, and the OpenMP multi-platform, shared-memory, parallel programing method is used to achieve real-time performance.