GPU-Accelerated Path-Dependent Marginal Information Gain for Autonomous Exploration

Consider an initially unknown static three-dimensional environment $\mathcal{W}\subset\mathbb{R}^{3}$ and a UAV equipped with a depth camera. The goal of an exploration planner is to generate collision-free feasible paths that build a complete map of $\mathcal{W}$ from the measurements acquired by the perception sensor. Since the environment is not known a priori, exploration is performed online in a receding-horizon fashion, replanning as new observations become available.

At each planning iteration, the current map $\mathcal{W}_{\mathrm{map}}$ partitions the environment into voxels of edge length $v_{\mathrm{size}}$, each labelled as free, occupied, or unknown. A sampling-based planner uses $\mathcal{W}_{\mathrm{map}}$ to construct a tree $\mathcal{T}=(\mathcal{N},\mathcal{E})$ rooted at the UAV’s current state, where each edge $e\in\mathcal{E}$ is a collision-free segment connecting two nodes. In turn, each node $n\in\mathcal{N}$ is characterized by a candidate viewpoint position $\mathbf{p}_{n}\in\mathbb{R}^{3}$ and yaw $\psi_{n}$, information gain $g(n)$, cumulative traversal cost $c(n)$, and path objective $o(n)$. The planner ranks the candidate paths according to $o(n)$, executes the first segment of the highest-valued path, and replans.

For a candidate position $\mathbf{p}_{n}$ and yaw $\psi$, let $\mathcal{I}(\mathbf{p}_{n},\psi)$ be the set of unknown voxels visible from that viewpoint. Existing approaches [3, 4, 6, 5] evaluate candidate viewpoints independently using absolute gain, defined as

where $v_{\mathrm{size}}^{3}$ is the voxel volume and $\left|\mathcal{I}(\mathbf{p}_{n},\psi)\right|$ denotes the number of visible unknown voxels.

The yaw $\psi_{n}$ is selected by maximizing $g_{\mathrm{abs}}(\mathbf{p}_{n},\psi)$ over $\psi\in\Psi$, and the resulting gain is stored as $g_{\mathrm{abs}}(n)$. Since this gain depends only on the current viewpoint, overlap with preceding viewpoints along the same path is ignored. The same unknown voxels may contribute to multiple nodes, overestimating the accumulated path gain. As overlap differs between paths, this can also change their relative ranking.

To eliminate the mutual independence assumption, the gain of a node must be conditioned on the observations expected from the previous viewpoints along its path. Let $\mathrm{Pa}(n)$ be the parent of node $n$ and $\mathcal{A}(n)=\{\mathrm{Pa}(n),\mathrm{Pa}(\mathrm{Pa}(n)),\ldots,n_{\mathrm{root}}\}$ its ancestor set, containing all preceding nodes along the path to $n$. Marginal gain accounts for all ancestors,

Similar to the absolute gain, the yaw $\psi_{n}$ is selected by maximizing $g_{\mathrm{all}}(\mathbf{p}_{n},\psi)$ over $\psi\in\Psi$, and the resulting gain is stored as $g_{\mathrm{all}}(n)$. Both the gain and yaw depend on the ancestor observations and selected yaws. The nodes must therefore be evaluated in tree-depth order. The two formulations satisfy $g_{\mathrm{all}}(n)\leq g_{\mathrm{abs}}(n)$.

Candidate paths are ranked by an objective function that combines information gain and traversal cost. The cumulative cost $c(n)$ to node $n$ is

corresponding to the traveled distance along the path. We use the exponentially decaying objective commonly adopted in sampling-based planners [3, 4, 5],

where $\lambda>0$ is a penalizing coefficient for traversal cost. This favors informative paths while reducing the contribution of viewpoints that are farther from the tree root.

Each node $n$ stores its viewpoint observation as a depth buffer $D_{n}\in\mathbb{R}^{W_{D}\times H_{D}}$, replacing the voxel hash map used by an exact implementation. Each pixel stores the camera-frame $z$-depth of the first occupied voxel along its ray, or of the ray endpoint at the maximum sensing range $d_{\max}$ if no occupied voxel is hit. Free and unknown voxels are traversed without terminating the ray.

The buffer resolution is chosen so that one pixel at the maximum sensing range matches the voxel size $v_{\mathrm{size}}$. For horizontal and vertical fields of view $\theta_{h}$ and $\theta_{v}$, we set $W_{D}=\left\lceil 2d_{\max}\tan(\theta_{h}/2)/v_{\mathrm{size}}\right\rceil$ and $H_{D}=\left\lceil 2d_{\max}\tan(\theta_{v}/2)/v_{\mathrm{size}}\right\rceil$. Once $\psi_{n}$ is selected, its depth buffer is kept on the GPU and used to evaluate descendants.

Information gain is evaluated by raycasting from $(\mathbf{p}_{n},\psi)$ in the current map. For marginal gain, the portions of each ray already observed by the ancestors of $n$ are identified using their depth buffers. To determine this overlap, the candidate ray is projected into each ancestor buffer and its depth is compared with the stored values. The camera is modeled using a pinhole camera model to perform this comparison in screen space.

Consider a candidate ray $\mathbf{r}(t)=\mathbf{p}_{n}+t\mathbf{d}$, with $0\leq t\leq d_{\max}$, where $\mathbf{d}$ is its direction. For each ancestor $a\in\mathcal{A}(n)$, the ray origin and direction are transformed into the ancestor camera frame and the ray is clipped along the camera $z$-axis to the valid depth range. The resulting 3D segment is projected onto the ancestor depth buffer as

where $f_{x}$ and $f_{y}$ are the focal lengths and $(c_{x},c_{y})$ is the principal point. The resulting 2D segment is clipped to the depth buffer bounds using the Liang-Barsky algorithm [17]. Together, the depth and image clipping restrict the overlap check to the valid camera frustum and, with the buffer resolution defined previously, keep the spacing between samples of the projected ray smaller than or equal to voxel size. Without this clipping, segments outside the sensing range or field of view would be undersampled in screen space, introducing errors in the depth interpolation and overlap estimation.

The clipped segment is traversed using a 2D digital differential analyzer (DDA) [8]. Under perspective projection, depth is nonlinear in screen space, but inverse depth $w=1/z$ is linear. So, for interpolation parameter $s\in[0,1]$,

where $w_{0}$=$1/z_{0}$ and $w_{1}$=$1/z_{1}$ are the inverse depths at segment endpoints. At each visited pixel, $z(s)$ is compared with the depth in $D_{a}$, and the segment is considered previously observed when $z(s)$ is smaller or equal to the stored depth.

Since this condition can change several times along the segment, one ancestor can produce multiple overlap intervals,

where $K_{a}$ is the number of overlap intervals from ancestor $a$. The ancestors are processed sequentially, merging each new $\mathcal{S}_{a}(\mathbf{r})$ with the accumulated intervals when they overlap.

Once the overlap intervals are determined, a 3D DDA computes the gain. Voxels inside the merged intervals are traversed but do not contribute, while unknown voxels outside them contribute normally and occupied voxels terminate the ray traversal. Figure 3 illustrates this process.

Marginal gain introduces an ordering constraint because the gain and selected yaw of a node depend on its ancestors. However, this constraint is limited to tree depth, since nodes at the same depth cannot be ancestors of one another. Consequently, once all shallower levels are processed, the nodes at the current depth can be evaluated independently.

Let $\mathcal{L}_{d}=\{n\in\mathcal{N}:\operatorname{depth}(n)=d\}$ be the set of nodes at depth $d$. The GPU evaluates the tree in increasing depth order. For each $n\in\mathcal{L}_{d}$, all candidate yaws are evaluated against the depth buffers of its ancestors and the yaw maximizing the gain is selected. Nodes and sensor rays within the same level are evaluated in parallel and, after each level, the resulting gains and yaws are returned to the CPU.

Therefore, the evaluation is sequential across depth levels but parallel within each level, preserving the ancestor ordering of marginal gain while still exploiting GPU parallelism.

Conventional sampling-based planners evaluate each node as it is generated. However, this is not suitable for a batched, depth-ordered evaluation like the one proposed. Therefore, in our framework, the CPU first builds the complete candidate tree by sampling and collision checking, without evaluating gain during tree expansion. The tree is then organized by depth and evaluated on the GPU, which returns the gain and selected yaw of each node to the CPU. Finally, the CPU computes the traversal costs and path objectives and selects the best branch to execute. By delaying gain evaluation until the tree is complete, enough candidates are provided to use the GPU parallelism. This structure is integrated into RH-NBVP [3] and AEP [4], whose sampling, collision checking, path scoring, branch selection and replanning logic are unchanged. Only the per-node gain evaluation is replaced by the depth-ordered GPU evaluation of the complete tree.

The 3D DDA used to compute information gain is common to both the exact voxel hash map implementation and the proposed depth buffer method. As such, we only compare the additional time and space complexity needed to condition a candidate observation on its ancestors. Let $A_{n}=|\mathcal{A}(n)|$ be the number of ancestors of node $n$. In the hash map implementation, for a fixed observed volume, the number of unknown voxels in an observation scales as $\mathcal{O}(v_{\mathrm{size}}^{-3})$. Evaluating a candidate requires merging the $A_{n}$ ancestor maps and checking its voxels against the merged set. Assuming $\mathcal{O}(1)$ average insertion and lookup, the time complexity is $T_{\mathrm{hash}}=\mathcal{O}(A_{n}v_{\mathrm{size}}^{-3})$. Each stored observation requires $\mathcal{O}(v_{\mathrm{size}}^{-3})$ space, while the merged ancestor map can grow to $\mathcal{O}(A_{n}v_{\mathrm{size}}^{-3})$. To avoid storing a merged map per node, only local observations are stored, and the merged map is rebuilt when needed.

The proposed method stores one depth buffer per node. With $N_{r}=W_{D}H_{D}$ rays and $W_{D},H_{D}=\mathcal{O}(v_{\mathrm{size}}^{-1})$, each observation requires $N_{r}=\mathcal{O}(v_{\mathrm{size}}^{-2})$ depth values. Each candidate ray is checked against the $A_{n}$ ancestor buffers and traverses at most $L_{D}=\mathcal{O}(W_{D}+H_{D})=\mathcal{O}(v_{\mathrm{size}}^{-1})$ pixels in each buffer. The time complexity is $T_{\mathrm{depth}}=\mathcal{O}(A_{n}N_{r}L_{D})=\mathcal{O}(A_{n}v_{\mathrm{size}}^{-3})$.

Both approaches have the same worst-case time complexity, but differ by a factor of $\mathcal{O}(v_{\mathrm{size}}^{-1})$ in space complexity, with the hash map additionally constructing merged ancestor maps during evaluation. A GPU implementation of the hash map approach is possible, but the path-dependent formulation requires a merged set for each candidate with up to $\mathcal{O}(A_{n}v_{\mathrm{size}}^{-3})$ entries, whose size is unknown in advance and which must be written by multiple threads at once. Our method avoids this, since its $\mathcal{O}(v_{\mathrm{size}}^{-2})$ candidate rays are independent and can be evaluated in parallel, while the $\mathcal{O}(A_{n}v_{\mathrm{size}}^{-1})$ traversal stays local to each ray.

We first validate the proposed GPU depth buffer marginal gain against the exact CPU hash map implementation using fixed yaw angles so that both evaluate the same viewpoints. Across 166,440 node evaluations, absolute gain matches the reference exactly, while the proposed marginal gain achieves $R^{2}=0.9937$ with a regression slope of $1.028$. The resulting difference comes because each depth buffer pixel represents the area covered by a beam instead of a ray, and can correspond to more than one voxel, up to four at the resolutions used. When the depth is interpolated, this approximation error accumulates along the beam, resulting in a small overestimation of the gain. Figure 4 shows that the proposed method remains close to the exact marginal gain, with a $5$-$10\%$ overestimation. In contrast, single-parent $g_{sp}$ and absolute gain errors increase with branch depth as more ancestor observations accumulate, reaching up to $2\times$ and $2.6\times$ the exact marginal gain, respectively.

Table II compares the CPU and GPU evaluation times of absolute and marginal gain for different tree sizes $N$ on the desktop and Jetson Orin NX. GPU timings include gain computation and CPU-GPU data transfers and are averaged over 10 replanning iterations. The Jetson timings were obtained using hardware-in-the-loop, with Gazebo running on the desktop while mapping, planning, and gain evaluation ran on the Jetson. On the desktop, the proposed marginal gain is $6$-$48\times$ faster than the hash map implementation at $0.2\text{\,}\mathrm{m}$ and $12$-$118\times$ faster at $0.1\text{\,}\mathrm{m}$. For $N=10000$, evaluation time drops from $16.30\text{\,}\mathrm{s}$ to $0.40\text{\,}\mathrm{s}$ and from $181.57\text{\,}\mathrm{s}$ to $1.95\text{\,}\mathrm{s}$, respectively. The Jetson shows the same trend, achieving speed-ups of $3$-$12\times$ and $10$-$28\times$. The method also scales better with map resolution since halving the voxel size increases the geometric mean marginal-gain time by $14.3\times$ on the desktop CPU versus $5.2\times$ on the GPU, and by $10.7\times$ versus $4.4\times$ on the Jetson. GPU absolute gain remains faster than marginal gain in every configuration, so any exploration improvement reported in the following sections is obtained despite a higher gain evaluation cost.

We compare AEP and RH-NBVP in the school environment, varying $N_{\max}$ from 50 to 500, with $N_{\mathrm{termination}}$ scaled accordingly. Table III summarizes the results. For both planners, marginal gain improves exploration efficiency, reducing $E_{95}$. With RH-NBVP, the effect becomes more significant as the tree grows. At $N_{\max}=50$ the difference in $E_{95}$ is small, since the limited number of branches often leads both gains to select the same direction. Larger trees provide more alternative branches, making overlap more likely to affect their ranking and change the selected branch. AEP benefits from marginal gain even with the smallest tree. Removing overlapping regions makes low-information branches fall below $g_{\mathrm{zero}}$ earlier, allowing AEP to switch to its global planner sooner. It also improves the decisions of the global planner, which, unlike RH-NBVP that only executes the first segment of the selected branch, executes the full path to the selected waypoint. Therefore, accounting for overlap along that path improves the waypoint selection. This is also reflected in the shorter paths and earlier termination times obtained across all tree sizes.

We compare absolute and marginal gain using RH-NBVP and AEP across the three environments, with the exploration curves in Figure 5. AEP reaches every coverage milestone earlier with marginal gain in all environments, and the improvement grows with coverage. Between $E_{25}$ and $E_{95}$ it rises from 0.8 to 1.8 min in the school, from 1.3 to 2.3 min in the large maze, and from 0.2 to 1.1 min in the multi-story environment. This is not explained by faster motion, since the average velocities are similar for both formulations. Marginal gain completes exploration with a shorter path in two of the three environments, from 895 to 799 m in the school and 520 to 505 m in the multi-story environment. The large maze keeps a similar traveled distance, at 679 against 684 m, while $E_{95}$ decreases by 2.3 min. AEP also terminates earlier in every environment, by 2.0 min in the school, and 0.5 min in the other environments. The improvement in $E_{95}$ exceeds the improvement in termination time in the large maze, indicating that most of it occurs before high coverage is reached. Afterward, only small and scattered regions remain, resulting in similar termination times.

For RH-NBVP, the improvement depends more on the environment, reducing $E_{95}$ by 0.7 min in the school, and 3.8 min in the large maze. The large maze shows the most consistent improvement, growing from 2.4 min at $E_{25}$ to 3.8 min at $E_{95}$. This is also visible in Fig. 6, where marginal gain increases the explored volume after $15$ min from $58.8\%$ to $75.0\%$. In the school the difference increases to 2.1 min at $E_{75}$ before decreasing to 0.7 min at $E_{95}$. The path length and velocity show no consistent advantage for either formulation.

The multi-story environment is the only case where marginal gain does not improve $E_{95}$ with RH-NBVP. Although it reaches $E_{25}$ and $E_{50}$ earlier, absolute gain becomes 0.4 min faster at $E_{75}$ and $E_{95}$. The tighter layout and the required $d_{\mathrm{collision}}$ make it harder for the sampling tree to reach informative regions, so additional nodes do not necessarily improve access to unexplored areas. The gain formulation can change the rank of reachable viewpoints, but it cannot compensate for limitations in the sampling process. The same effect is visible with AEP, where the termination time differs by only 0.5 min. Final coverage remains similar for both formulations in all environments, indicating that marginal gain improves exploration efficiency but not the final mapped volume.

We evaluate RH-NBVP in the large maze with execution horizons $h=1$, $3$, and $5$, where $h$ is the number of consecutive nodes executed before replanning. The large maze is used because its size requires many planning decisions throughout the mission, making the effect of committing to longer branches more evident. Table IV reports the time $E_{95}$, path length $L_{95}$ and average velocity $\bar{v}_{95}$ to reach $95\%$ coverage, together with the total planning time $T_{\mathrm{plan}}$.

The results show that $E_{95}$ remains similar as the horizon increases with marginal gain. Longer horizons require fewer replanning iterations, so $T_{\mathrm{plan}}$ decreases from $353\text{\,}\mathrm{s}$ at $h=1$ to $105\text{\,}\mathrm{s}$ at $h=5$, while $\bar{v}_{95}$ increases from $0.35\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ to $0.42\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$. However, path length increases with the execution horizon. With $h=1$, the planner updates its decision after every executed node using the latest map, resulting in more efficient paths. For longer horizons, the UAV must execute more of the selected branch before changing direction. As a result, $L_{95}$ with marginal gain grows from $565\text{\,}\mathrm{m}$ to $635\text{\,}\mathrm{m}$ and $659\text{\,}\mathrm{m}$, and absolute gain follows the same trend, from $635\text{\,}\mathrm{m}$ to $696\text{\,}\mathrm{m}$ and $743\text{\,}\mathrm{m}$. Therefore, longer horizons trade path efficiency for reduced replanning time and higher velocities, keeping $E_{95}$ similar.

More importantly, the benefit of marginal gain persists at longer horizons. At $h=3$ it reaches $95\%$ coverage with $635\text{\,}\mathrm{m}$, the same mean path length as absolute gain at $h=1$, needing just $24\text{\,}\mathrm{m}$ more at $h=5$. As such, marginal gain at $h=3$ and $h=5$ achieves similar path efficiency to absolute gain replanning after every node, also reducing $E_{95}$ from $30.5\text{\,}\mathrm{min}$ to $26.0\text{\,}\mathrm{min}$ and $26.5\text{\,}\mathrm{min}$. At $h=3$, this is achieved with comparable planning times, $151\text{\,}\mathrm{s}$ versus $156\text{\,}\mathrm{s}$.