From Topological Distance to Spatial Buffers: Efficient Halo Queries on Triangular DGGS

1. Introduction

Discrete global grid systems (DGGSs) provide a geospatial platform [1,2]. Systems such as H3 [3], S2 [4], and rHEALPix [5], together with research systems [6,7,8], partition the Earth into hierarchies of indexed cells at multiple resolutions. These cell indices support aggregation, joins, and neighborhood queries across applications such as urban analytics, Earth observation, climate modeling, and logistics.

Spatial proximity is an important operation in these applications. Buffering a road network, identifying facilities within a service radius, and retrieving parcels near a flood zone ask a common question: given a source region, what lies within a specified distance? In a DGGS, proximity also has a discrete form. Starting from cells covering the source, one expands through cells that share an edge and counts the adjacency steps. This count is the topological distance. Cell indices make these steps suitable for parallel execution because adjacent cells can be constructed without repeatedly measuring geographic distances. Two practical questions follow: how can topological rings be computed efficiently at fine resolutions, and how can their expansion support complete coverage of a requested distance in meters?

Our computational workload begins with the stored cell indices needed by the first question. A halo query constructs the cells reached after one, two, and up to a requested number of edge adjacency steps. Such queries support repeated neighborhood access to data already assigned to a DGGS, for example when changing the neighborhood depth around a fixed indexed feature. Their cost is distinct from assigning a new geographic dataset to cells or evaluating a buffer in meters. The measured acceleration concerns this indexed halo operation. The second question is addressed by a separate boundary buffer study that adds source covering and geometric distance tests around the same halo component.

A concrete example connects these questions. Consider all locations within 1 km of the complete boundary of Nose Hill Park in Calgary. The source is the boundary curve, not the filled park, so the requested region extends to both sides of that curve, as shown in Figure 1(a). Cells covering the complete curve form the seed set P. Starting from these cells, the halo query repeatedly finds cells that share an edge with the current ring. Cells first reached after one step form the halo ring$H_1$; cells first reached after two steps form $H_2$, and so on, up to the requested hop depth K. Figure 1(b) shows the first three rings, and Section 3.1 gives the formal definition. The same boundary request is used in the measured ellipsoidal experiment.

A ring count alone does not establish that the whole 1 km boundary buffer has been reached. The procedure must first cover the complete source, expand enough rings, and then exclude cells that cannot intersect the requested region. Figure 1(c) illustrates a conservative cover: every required cell is retained, but some additional cells may remain. Geometric subdivision can tighten a distance test for an uncertain cell without changing the cell being returned or the query resolution. Choosing a finer query resolution is a different operation that produces smaller returned cells. The coverage conditions and the numerical limits of filtering are developed in Section 3.5.1 and Section 3.5.2.

The distance request therefore has two connected stages: construct the cells reached by adjacency, then use geometric tests to decide whether expansion is sufficient and which cells to retain. Efficient ring construction is the shared computational foundation. The first computational difficulty is the size of the refined grid. A triangular base mesh, meaning the coarse polyhedron from which cells are refined, contains $F{4}^r$ cells at level r when each of its F faces is repeatedly divided into four child triangles. Fine resolutions therefore produce far more cells than a local query needs. Storing three adjacent cell indices for every refined cell would allocate memory for the whole Earth even when the requested halo covers only a small region. We instead need to construct adjacency only for cells involved in the query.

The second difficulty concerns how each cell is indexed. A hierarchical index records the base mesh face and one child digit for every refinement. This form supports direct access to ancestors and contiguous ranges of descendants [8,10]. Constructing an adjacent cell at the same resolution can instead require decoding and updating the geometric state represented by those digits. Integer coordinates support adjacency through a small number of arithmetic operations, but do not express the refinement sequence as a visible prefix. A halo query therefore benefits from a hierarchical index for storage and hierarchy operations, and a coordinate representation for repeated adjacency calculations.

Parallel expansion introduces the third difficulty. Each triangular cell in the current ring has three edge adjacent cells. Different workers can generate the same adjacent cell, and a generated cell may also belong to the seed set or an earlier ring. The query must reject these repeated discoveries and add every new cell exactly once. The central problem is therefore to construct exact halo rings at fine resolutions while storing query sized state and coordinating parallel discoveries, without storing the complete refined adjacency structure.

We address these three difficulties with connected design choices. First, a dual index architecture uses one index for the query interface and another during expansion. Cells enter and leave the query as a public hierarchical index, which records the base mesh face, refinement sequence, and resolution. Each input is converted once into a compact integer triangular index, which stores the face, two integer coordinates, and the triangle orientation in one 64 bit word. This representation constructs adjacent cells at the shared query resolution through direct arithmetic. Only the completed output is converted back. Section 3.3.3 places this representation in the context of existing coordinate indexes; the performance evaluation compares our execution variants on the same topology.

Second, adjacency is reconstructed rather than stored. Each base mesh face has a chart: a planar coordinate domain used for its triangular subdivision. Fixed integer rules construct adjacent cells that remain in the same chart. For a cell beside a base mesh edge, the Atlas of Connectivity Maps (ACM) records how coordinates continue into the neighboring chart [11]. Each record identifies the destination face, the shared edge, and the ordering of cells along that edge. Deterministic refinement allows the same records to apply at every resolution. Local arithmetic and $3F$ directed base mesh edge records therefore describe the refined adjacency without storing a record for each refined cell.

Third, the expansion state remains in graphics processing unit (GPU) memory. This state includes the seed cells, cells already discovered, and the current and next rings. A hash table stores the discovered cell indices. An atomic insertion allows only one worker to register a new index, even if other workers attempt the same insertion simultaneously. The successful worker adds the cell to the next ring; workers that find the index already present do not add it again. This connects duplicate removal directly to ring construction and avoids transferring intermediate rings to the central processing unit (CPU). Figure 5 shows the complete query sequence, while Figure 10 explains the insertion procedure after its data structures have been introduced.

The Pentakis dodecahedron is our primary research base mesh, following the triangular mesh construction studied by Wecker et al. [8]. We choose it to test the atlas on multiple face arrangements, not because of a demonstrated adoption or performance advantage over established DGGS libraries. Its 60 triangular faces include vertices where five faces meet and others where six meet. Testing both arrangements checks movement among different numbers of adjacent faces. At level 28, this mesh contains approximately $4.323\times {10}^{18}$ cells, with a global average area of about $1.18$ square centimeters per cell. The architecture nevertheless stores only its 180 directed base mesh edge records and allocated query state. Within the current six bit face field, the same execution method supports oriented triangular base meshes with at most 64 faces that use the same subdivision. An icosahedron requires 60 directed records and an octahedron requires 24. Their performance must still be evaluated separately under comparable spatial workloads.

We evaluate the architecture by comparing parallel performance, isolating component contributions, measuring complete workflow costs, and testing correctness. Official City of Calgary building, street, and boundary data supply 81 configurations through level 28. On the tested Intel Core i7 12700K and NVIDIA RTX 4070, the GPU implementation is faster than the better of two CPU baselines that store the discovered set in a flat open addressing table (the flat table baselines of Section 3.4.1) in 48 of the 81 configurations; the median ratio over all configurations is $2.38$, and the 33 configurations won by the CPU return at most 3,000 cells. The largest query constructs 10.8 million halo cells from 100,000 building seed cells in 13.3 ms, measured to page locked (pinned) host memory, compared with 117 ms for the twenty thread flat table. Sections 3.6.4 and 4.1.1 state the timing boundaries, the selection of measurement sessions, and its sensitivity. These ratios compare GPU pinned host readiness with CPU returned container completion; the additional GPU container copy is reported separately. Mapping points to cells, source covering, and rendering are excluded from these times. Independent topology tests exercise every directed base mesh edge and every base mesh vertex of the Pentakis mesh; output agreement checks for the measured configurations are reported separately.

The second stage uses these rings to evaluate the same boundary distance request. It constructs World Geodetic System 1984 (WGS84) buffers using complete source covering, adaptive halo expansion, and curved cell distance tests. In all seven tested buffer configurations, the CPU only variant is faster than the hybrid variant, which runs ring expansion on the GPU and the geometric tests on the CPU, because geometric processing and GPU initialization dominate the total cost. This is distinct from the measured acceleration of the reusable halo component. The implementation retains unresolved cells and states its numerical allowances. No formal correctness proof of the floating point distance implementation is provided; its implementation evidence is numerical validation.

The primary contribution is an execution architecture for exact halo construction from stored indices. The boundary buffer study demonstrates the additional geometric requirements and their present cost. The first stage addresses the cost of large adjacency expansions; the second specifies the additional geometric conditions needed for a distance request. The contributions are:

• Query specific representation and execution. A public hierarchical index serves storage and hierarchy operations; a compact integer triangular index serves repeated adjacency and discovery. Building on established indexing and atlas principles, the architecture reconstructs neighbors from arithmetic and $3F$ base mesh edge records, without refined connectivity arrays.
• Concurrent exact ring construction. A GPU resident discovered set rejects repeated cells and authorizes exactly one append for every newly admitted cell. Cooperative ring scheduling retains intermediate state on the device. Memory depends on input, hash, and output capacities rather than the global cell count.
• Independent validation. Geometric references that do not consult the runtime atlas check base mesh and refined adjacency. Globally distributed queries exercise all directed base mesh edges and base mesh vertices of the Pentakis mesh, while failure tests require explicit errors instead of incomplete results.
• Measured performance and limits. Serial and multithreaded flat table baselines, a thread count sweep, stage measurements, and controlled mechanism changes identify the workload and delivery conditions that support acceleration. Sensitivity to the measurement session and a local routing check state the limits of those findings.
• Distance coverage conditions. Complete source covering, a whole ring stopping argument, and cell extent bounds connect exact topological expansion to conservative spatial coverage. A separate WGS84 boundary study evaluates the additional geometric work and its numerical limits; CPU execution is faster for its tested complete queries.

The work does not introduce a new projection or replace established DGGS indexing principles. It connects representations suited to repeated adjacency with geometric tests suited to distance requests. After defining exact rings and their spherical interpretation, we explain the indices and execution methods used to construct them. We then add source covering, stopping, and filtering to obtain boundary buffers. The evaluation separates the performance of the halo component from the cost and numerical limitations of the complete query.

2. Background and Related Work

This section introduces the foundations needed to understand the problem and the proposed solution. We first review how a discrete global grid system divides the Earth, identifies its cells, and defines neighborhood operations. We then explain the triangular coordinates and connectivity atlases that make local cell adjacency possible. Finally, we review parallel frontier processing and identify the remaining difficulty: constructing a very large set of newly discovered cells without storing the complete refined topology.

2.1. DGGS Foundations: Domains, Indexing, and Neighborhood Operations

A discrete global grid system (DGGS) divides a spherical or ellipsoidal Earth model into cells at several resolutions. Its projection, refinement, and indexing rules determine how locations are assigned to cells and how those cells are organized [1]. The cell shape and the coarse polyhedron, called the base mesh, strongly influence refinement, distortion, and neighborhood structure. Hexagons provide nearly uniform adjacency, but their fractional refinement does not form a simple congruent subdivision. Quadrilaterals align naturally with raster storage; their distortion depends on the base mesh and projection, with polar distortion arising in some constructions. Triangles are planar primitives that admit a congruent subdivision in which every triangle is replaced by four smaller triangles [6,12]. Prominent triangular base meshes include the octahedron used by the Hierarchical Triangular Mesh [13], the icosahedron used by many geodesic grids [12], and the Disdyakis triacontahedron with 120 triangular faces [7]. We use the Pentakis dodecahedron with 60 triangular faces as the primary example. Each face supplies a planar coordinate domain, called a chart, for projection and subdivision. The smaller angular extent of these faces limits the extent of each projection chart. Five triangular faces meet at some vertices and six meet at others, allowing boundary transitions to be tested in both arrangements [8].

Once the base mesh and refinement define the cells, an index is needed to locate each cell and express its relationships to others. The indexing method determines which DGGS operations can be performed directly. A hierarchical prefix index records one refinement digit for every parent to child step. It therefore exposes the resolution and ancestry of a cell, but finding a same level adjacent cell can require decoding and updating geometric state at several refinement steps [10,14]. A space filling curve, such as Morton order [15] or Hilbert order, maps cells to a one dimensional sequence with useful locality, but cells that share an edge are not necessarily consecutive on the curve [10,16]. Conversely, a coordinate based index resolves same level cell adjacency in constant time through simple arithmetic, but obtains hierarchical relationships through operations such as integer division or bit shifting rather than through an explicit prefix. Constant time adjacency is not exclusive to coordinate indexes: Lee and Samet construct equal size neighbors in linear triangle quadtrees using a fixed number of bit operations [17]. Our distinction is therefore the representation used throughout a bulk query, not a claim that every hierarchical neighbor method is slower. The dual index architecture uses both strengths. The query interface accepts and returns a public hierarchical index for persistent storage and hierarchy operations. At the beginning of a query, each public hierarchical index is converted once into a compact integer triangular index. All repeated cell adjacency operations use this compact form, and only the final result is converted back.

These indexing choices matter because spatial operations repeatedly access cells and their neighborhoods. DGGS operations include multiresolution data integration, data retrieval, spatial predicates, aggregation, buffering, and distance analysis [1,2,18]. Kazemi et al. use a coarse to fine distance transform to limit the region evaluated at each refinement level [19]. Alderson et al. distinguish definitely included and uncertain cells when offsetting spherical curves, and refine boundary cells to improve precision [20]. Lansang and Samavati explicitly distinguish discrete sample coverage from continuous coverage, using a bound on cell extent to restore the latter [21]. These studies show why adjacency count and a spatial buffer are related but different operations. Our halo query constructs exact discrete rings; the additional distance formulation conservatively covers a specified spatial region. At high resolution, neither operation should require storage for all $F{4}^r$ refined cells and their adjacency records when only a smaller region is requested.

The evaluation isolates the execution method by requiring serial central processing unit (CPU), multithreaded CPU, and graphics processing unit (GPU) implementations to process the same Pentakis cells with identical adjacency and output definitions. rHEALPix, H3, and S2 are important DGGS designs with different cell shapes, refinements, and neighborhood degrees [3,4,5]. Related library operations include H3 grid disks and rings and S2 cell neighbors [22,23]. Their cell degrees and single cell or set semantics must be aligned before timing comparisons are meaningful. Zheng et al. instead accelerate raster assignment to a hexagonal DGGS on the GPU [24]; that preparation operation differs from expansion of already indexed cells. Volumetric DGGSs add a radial index and vertical adjacency to the surface index [25]. A comparison among these systems would mix execution cost with differences in geometry and neighborhood semantics. Such a comparison is outside the present experimental scope. The experiments instead hold the DGGS topology fixed so that the measured differences describe the execution architecture.

2.2. Triangular Indexing and Connectivity Atlases

The previous subsection identifies same level cell adjacency as the operation repeated most often in a halo query. We now review the triangular coordinate systems and boundary transformations that make this operation possible. The Barycentric Indexing Method (BIM) assigns a symmetric three axis integer coordinate $(a,b,c)$ to elements of a refined base triangle. A complete BIM index also records the base cell and resolution, and supports hierarchy, neighborhood access, and movement across triangular meshes [8]. Earlier work on Hexagonal Connectivity Maps showed that coordinates aligned with a diamond hierarchy can support constant time neighborhood and hierarchy operations on refined triangular polyhedra [26]. These methods establish an important principle: the internal coordinate representation should expose the operations that dominate the intended workload.

Local coordinates describe movement inside a face, but movement across a face boundary also requires a relation between the two coordinate systems. A Connectivity Map records how regularly subdivided cells connect inside one chart. An Atlas of Connectivity Maps (ACM) extends this idea to several charts by storing the finite coordinate transformations required when a cell adjacency crosses a base mesh edge [11,27,28]. Related work applies these structures to adaptive connectivity and to the storage, retrieval, transmission, and compression of multiresolution geospatial fields [18,28]. The literature therefore supplies two foundations used in this paper: integer coordinates describe regular adjacency inside a triangular chart, and an atlas reconciles those coordinates at chart boundaries. Previous studies primarily examine individual hierarchy or neighborhood operations and multiresolution field organization. Our work extends this view to a complete halo frontier generated simultaneously from many seed cells.

2.3. GPU Frontier Processing

Integer coordinates and an atlas make the three adjacent cells of one triangular cell inexpensive to obtain. A halo query, however, expands many cells simultaneously, and two expanded cells can discover the same adjacent cell. This subsection reviews the parallel methods needed to form the next ring without duplicates. In each iteration, the frontier is the completed halo ring being expanded. Each frontier cell generates three candidate cells, one for each shared edge. Work efficient GPU frontier methods organize active elements into such frontiers and suppress repeated visits [29,30]. General adjacency frameworks load stored neighbor lists. A triangular DGGS does not need those lists because the three candidates can be reconstructed from the cell index.

The method must remember which candidate cells have already been discovered because adjacent frontier cells can generate the same candidate. Spatial workloads commonly use a hash table to represent a sparse discovered set without allocating an entry for every possible cell [31]. Concurrent GPU hash tables show that groups of GPU threads and atomic updates can provide high insertion throughput [32]. GPU hashing is an established parallel primitive rather than a novelty of this work; cuCollections, for example, provides GPU data structures based on concurrent operations [33]. An atomic compare and swap operation changes an empty table entry only if it still contains the expected empty value. It therefore allows one thread to claim a new cell while other threads that generated the same cell observe that it has already been claimed.

Concurrent insertion addresses duplicate discovery, but the next ring must also remain available for the following iteration. For this broader processing sequence, WarpCore demonstrates that cooperative probing and GPU hash table operations can provide high insertion and retrieval throughput while keeping a larger processing sequence in GPU memory [34]. Our architecture uses a more restricted set with no deletion during a query. A successful first insertion means that the candidate is a new halo cell and can be appended to the next frontier. An existing seed cell is excluded, a cell found in an earlier ring remains excluded, and a second insertion of the same candidate in the current iteration is suppressed as a duplicate. This direct relationship between insertion and halo membership allows one atomic decision to replace several separate set operations.

3. Materials and Methods

A complete boundary distance query follows four operations: cover the supplied curve with cells, construct complete topological rings, test whether expansion is sufficient, and filter the discovered cells by distance. The halo component supplies the second operation; source covering and geometric tests supply the others. This section follows their dependencies. Section 3.1 defines the two requested outputs, and Section 3.2 establishes the spherical domain and spatial interpretation. Section 3.3 explains the index representations and adjacency rules, and Section 3.4 explains their use in serial, multithreaded, and GPU ring construction. Section 3.5 then places source covering, stopping, and filtering around that component. Finally, Section 3.6 defines the studies used to evaluate each stage and the complete boundary query.

3.1. Topological Query Model and Problem Definition

Before choosing how to represent or process cells, we define the result that a halo query must return. We assume a closed triangular surface without boundary edges. Its base mesh is consistently oriented, meaning that adjacent faces use opposite directions along their shared edge. Each base mesh face is recursively divided into four smaller triangles, so every refined cell has exactly three edge adjacent cells. The following definition uses only these cells and their adjacency. The index representations and execution methods are introduced separately in Section 3.3 and Section 3.4.

The input feature is represented by a seed set P of cells at one resolution. The kth halo ring, denoted by $H_k$, contains the cells first reached after k adjacency steps. Here a ring is a set at equal minimum topological distance; it need not be connected or form one geometric annulus. This step count is topological distance. Spatial distance is instead measured in meters on the chosen Earth model; it does not determine membership in a fixed hop query. A preindexed halo query starts from stored cell indices, so geographic point assignment and source covering are separate preparation operations.

A distance query instead starts from a continuous source $\mathcal{S}$ and a distance d. It requests a cell cover of all locations within d of that source, so its required hop depth is not specified in advance. For the running example, $\mathcal{S}$ is the Nose Hill Park boundary. Section 3.5 defines its geometric target and coverage requirements after the cells and reusable halo component have been established.

3.1.1. Topological Query Model

We first define halo membership independently of geometry, indexing, and hardware. The seed set occupies topological distance zero; each later ring contains exactly the cells first reached after one additional edge adjacency step. The following notation makes that set semantics precise.

Let ${\mathcal{T}}_r=(C_r,E_r)$ denote the edge adjacency topology of a level r triangular DGGS, where $C_r$ is the complete set of cells at level r and

$$E_r=\left\{\{u,v\}\subset C_r:u\mathrm{and}v\mathrm{share}\mathrm{an}\mathrm{edge}\right\}$$

contains their edge adjacency relationships. Every refined cell is triangular and has three edge adjacent cells. Duplicate input indices are removed to form the seed set $P\subseteq C_r$, and the requested hop depth is a nonnegative integer K. For a set $S\subseteq C_r$, let $\Gamma \left(S\right)$ be the set of all cells that share an edge with at least one cell in S. Define the zeroth ring as $H_0=P$. The remaining disjoint rings follow the recurrence

$$H_k=\Gamma \left(H_{k-1}\right)\setminus {\displaystyle \bigcup _{q=0}^{k-1}}{H}_q,1\le k\le K.$$

(1)

Thus, for every $0\le k\le K$, $H_k$ contains exactly the cells whose minimum topological distance from P is k. The returned halo excludes the seed ring $H_0$ and is $H={\bigcup }_{k=1}^K{H}_k$, accompanied by ring offsets identifying the individual rings. When $K=0$, no adjacency expansion is requested and the returned halo is the empty set, while $H_0=P$ remains the seed ring. Duplicate input public hierarchical indices do not affect any ring because they are removed when constructing P.

3.1.2. Halo Construction from a Seed Set

The recurrence in Equation (1) defines which cells belong to each ring. We now express the same definition as a general algorithm before discussing any particular hardware. The seed set forms the ring $H_0$. At iteration k, the method collects every cell that shares an edge with a cell in the completed ring $H_{k-1}$. These collected cells are the candidates for the next ring. Candidates already present in $H_0,\dots ,H_{k-1}$ are removed. The remaining cells form the next disjoint ring $H_k$. The discovered set includes the seed set and every completed ring, while the returned halo excludes the seed set. This distinction ensures that each returned cell is assigned to its minimum topological distance.

Figure 2 builds the intuition before the formal procedure. The left panel expands the yellow seed ring $H_0=P$ and obtains the complete green ring $H_1$. The right panel expands that completed ring to obtain the light blue ring $H_2$, while the red cells highlight representative frontier members whose adjacent candidates contribute to the next set. The colors denote set level roles only; the parallel candidate classification mechanism used by the device implementation is introduced later in Section 3.4.4.

Algorithm 1 states this breadth first set expansion independently of any serial, multithreaded, or GPU realization.

Algorithm 1 Topological halo construction from a seed set

Require:   seed set $P\subseteq C_r$; hop depth $K\ge 0$ Ensure:   disjoint rings $H_1,\dots ,H_K$ and halo H 1: $H_0\leftarrow P$ 2: $V\leftarrow P$      ▹ all cells discovered so far 3: for $k\leftarrow 1$ to K do 4:     $Q\leftarrow \Gamma \left(H_{k-1}\right)$     ▹ adjacent candidates 5:     $H_k\leftarrow Q\setminus V$    ▹ new cells at distance k 6:     $V\leftarrow V\cup H_k$ 7: end for 8: return $H={\bigcup }_{k=1}^K{H}_k$ and the ring offsets

alg:halo expresses the query independently of any CPU or GPU realization. The discovered set maintains the invariant $V={\bigcup }_{q=0}^{k-1}{H}_q$ at the start of iteration k. The candidate set Q contains every cell adjacent to the previous ring, and subtracting V removes the seed cells and all cells discovered in earlier rings. Consequently, the remaining cells form exactly $H_k$. The update of V establishes the invariant for the next iteration. When $K=0$, the loop is skipped and the returned union is empty. Figure 2 visualizes the first two iterations at this set level: yellow denotes $H_0$, green denotes $H_1$, and light blue denotes $H_2$. The later GPU resident algorithm implements these same set transitions concurrently; its atomic candidate states are introduced separately in Section 3.4.4.

3.1.3. Formal Constraints and Design Goals

The recurrence specifies the correct output, but it does not yet express the scalability conditions that motivate the architecture. We now state those constraints formally so that each method component can be evaluated against a specific requirement.

Problem statement. Given the public hierarchical indices of P, their common resolution level r, and a nonnegative hop depth K, compute the exact disjoint sets $H_1,\dots ,H_K$, or the empty halo when $K=0$, without a refined global adjacency table. The static adjacency state is $\Theta \left(F\right)$ and contains exactly $3F$ directed base mesh edge transitions. Constructing one adjacent cell takes $O\left(1\right)$ time. If $n_{\mathrm{in}}$ is the supplied input count, M the allocated hash capacity, and B the allocated output capacity, dynamic storage is $O(n_{\mathrm{in}}+M+B+K)$. There is no required term proportional to $|C_r|$. This becomes $O\left(\right|P|+|H|+K)$ when input and chosen capacities are within constant factors of the realized sets; a pessimistic capacity estimate may allocate more. During accelerated execution, the state required to continue from one ring to the next remains in GPU memory until the complete result is materialized.

The design goals are:

• Algorithmic exactness. The architecture must evaluate the discrete expansion in Equation (1) exactly, without substituting approximate distance metrics, geometric envelopes, or incomplete candidate sets.
• Resolution scalability. The public hierarchical index must encode the level and refinement sequence and support direct ancestor access through level 28, while same level neighbor construction retains constant cost and introduces no storage term that grows with $|C_r|$.
• Base mesh portability. Boundary transitions must preserve exact topological adjacency for an oriented, watertight triangular base mesh that uses the same refinement rule. The current six bit face field supports meshes with at most 64 base mesh faces. Within this supported class, changing the polyhedron modifies its $\Theta \left(F\right)$ transition atlas but not the local expansion or parallel discovery procedure; a mesh with more faces requires a wider face field.
• Bounded failure. Invalid indices, insufficient hash capacity, excessive output, and unavailable memory must produce explicit errors rather than incomplete answers.

3.2. Spherical Geometry and Spatial Reach

Topological rings describe adjacency but do not yet give a distance in meters. We first define the spherical cells on which the query operates, then explain how centroid distances describe a completed ring. This geometric interpretation motivates the separate coverage tests in Section 3.5; the index representations used to perform the expansion follow in Section 3.3.

3.2.1. Base Mesh Geometry and Projection

The first part of the architecture defines the geometric domain needed for base mesh face ownership, chart orientation, and spatial measurements. The Pentakis dodecahedron has 32 vertices, 90 edges, and 60 triangular faces. Each face defines one triangular chart, meaning one planar coordinate domain associated with that face. At each refinement, midpoint subdivision replaces one triangle with four smaller triangles. Level r therefore has $N=2^r$ segments along each base mesh edge and $4^r$ planar subtriangles per face before radial normalization to the sphere. Five triangular faces meet at some Pentakis vertices and six meet at others. Testing both arrangements verifies boundary transitions among different numbers of adjacent base mesh faces. Figure 3 shows the primary base mesh; the two index representations are introduced separately in Sections 3.3.1 and 3.3.2.

Geographic point assignment uses a face wise gnomonic construction: a longitude and latitude position is mapped to the spherical domain, centrally projected onto a candidate base mesh face plane, tested by barycentric containment, and discretized into the refined triangular chart. Deterministic ownership resolves positions on shared base mesh edges. The central projection has an important line preserving property: a great circle plane through the sphere center intersects a base mesh face plane in a straight line. Consequently, straight refined edges in the planar chart map back under radial normalization to great circle geodesic arcs, and the boundaries of the resulting spherical triangular cells are geodesic segments. Cell centroids are likewise reconstructed on the face plane and radially normalized to the sphere. The gnomonic construction is neither distance preserving, equal area, nor conformal, but its distortion belongs to the underlying mesh based DGGS rather than to the halo architecture. More importantly for the present method, projection is an indexing stage operation and is absent from the repeated preindexed adjacency expansion shown in Figure 5.

Figure 3. The Pentakis dodecahedron used as the primary base mesh. Its 60 triangular faces define the planar charts from which refined spherical cells are constructed; two adjacent charts are highlighted in gold. The mesh has 32 vertices and 90 edges, corresponding to 180 directed base mesh edge transitions. Opaque face shading reveals the polyhedron’s three dimensional form.

Figure 3. The Pentakis dodecahedron used as the primary base mesh. Its 60 triangular faces define the planar charts from which refined spherical cells are constructed; two adjacent charts are highlighted in gold. The mesh has 32 vertices and 90 edges, corresponding to 180 directed base mesh edge transitions. Opaque face shading reveals the polyhedron’s three dimensional form.

3.2.2. Relating Topological Rings to Spatial Reach

A halo ring at topological distance k does not lie at a fixed spatial distance from the seed cells. Even on a uniform planar triangular grid, a path of k adjacency steps can turn in any direction, so topological distance k gives only an upper bound on spatial displacement, not a fixed multiple of a step length. On a DGGS the relationship is looser still, because projection distortion changes cell dimensions and shapes with resolution and with location across the base mesh faces. Exact distances between cells could be evaluated on the WGS84 ellipsoid with geodesic routines that account for the flattening of the Earth [35]. This subsection instead works on the spherical DGGS: cell centroids are radially reconstructed onto a sphere, and distances are evaluated as great circle arcs with one trigonometric operation per pair. The spatial extent of a completed ring $H_k$ is therefore characterized by the interval between the nearest and farthest such spherical centroid distances to the seed cells, as illustrated in Figure 4. These bounds are relative to the spherical model. Converting them to certified WGS84 intervals would require bounds that account for the adopted sphere, latitude, and direction; with the mean spherical radius used here and WGS84 flattening ($f\approx 1/298.257$), the local discrepancy near the poles is approximately $0.45\%$ relative to ellipsoidal distance. We retain the spherical bounds as a description of reach; whether K hops suffice to cover every location within a requested distance d requires both a coverage argument and explicit numerical limits, developed in Section 3.5.1.

The distances described above are between reconstructed cell centroids, not minimum distances between complete cell footprints. The quoted percentage describes the local polar limit for the adopted radius $R=6,371,008.8$ m; it is not a global error bound. At a pole, the WGS84 curvature radius is $a/(1-f_E)$, giving a relative local difference of $1-R(1-f_E)/a\approx 0.00447$, with a and $f_E$ defined in Section 3.5.2. For a finite route the discrepancy depends on its locations and directions. It is not used in any coverage certificate. Section 3.5.2 retains the same cell adjacency but evaluates distances directly on WGS84 and bounds the corresponding mapped cell footprints. Its separate ellipsoidal tests do not rely on converting spherical centroid distances by a fixed percentage.

Figure 4 illustrates the distinction for one pair of cells and one adjacency route. Repeating the spherical distance calculation over every cell in a completed ring gives the nearest and farthest distances from the seed set. These ring measurements describe outward reach; the coverage test in Section 3.5.1 must also rule out required cells that have not yet been reached.

Here $P\subseteq C_r$ is the canonical seed set: the level r DGGS cells representing the input feature after duplicate public hierarchical indices have been removed. A symbol $p\in P$ denotes one seed cell, whereas $c\in H_k$ denotes one cell in the kth halo ring. Let ${\mathbf{x}}_c,{\mathbf{x}}_p\in {\mathbb{S}}^2$ be the unit vectors of their radially reconstructed spherical centroids, and let R be the adopted spherical Earth radius. Their cell to cell spatial distance is the minor arc great circle geodesic length

$$d_{\mathbb{S}}(c,p)=2Ratan2\left(\parallel {\mathbf{x}}_c-{\mathbf{x}}_p\parallel ,\parallel {\mathbf{x}}_c+{\mathbf{x}}_p\parallel \right).$$

(2)

This form avoids subtracting nearly equal cosine values at centimeter separations. It is mathematically equivalent to the usual great circle formula, but a numerically stable formula alone is not a certified error bound. Using Equation (2), the distance from a ring centroid to its nearest seed centroid is

$$\delta (c,P)=\underset{p\in P}{\min }{d}_{\mathbb{S}}(c,p).$$

(3)

Here P is the same seed set defined in Section 3.1. This distance concerns the discrete seed centroids. It is not generally the distance to the original continuous source curve or to an entire seed cell.

For every nonempty topological ring $H_k$, the nearest and farthest ring centroids give the exact query specific spatial interval

$$L_k=\underset{c\in H_k}{\min }\delta (c,P),U_k=\underset{c\in H_k}{\max }\delta (c,P).$$

(4)

By Equation (4), every centroid in $H_k$ therefore lies between $L_k$ and $U_k$ meters from its nearest seed centroid. For a nonempty cumulative halo, its maximum centroid reach is ${\max }_{1\le k\le K,H_k\ne ⌀}{U}_k$. If the returned halo is empty, including when $K=0$, spatial reach is reported as not applicable rather than taking a maximum over an empty set. These reductions annotate a completed ring and do not alter halo membership.

A bound on outward reach is useful even when the exact ring intervals are not evaluated. Define

$${\ell }_r^{+}=\underset{\{u,v\}\in E_r}{\max }{d}_{\mathbb{S}}(u,v)$$

(5)

as the longest spatial step permitted by one legal adjacency move anywhere in the level r grid, measured between the spherical centroids of adjacent cells. Every $c\in H_k$ has an adjacency path of length k from at least one seed. The spherical triangle inequality therefore gives

$$0\le L_k\le U_k\le k{\ell }_r^{+}.$$

(6)

Together, Equations (5) and (6) bound where an already reached ring can lie. A sampled maximum is not a certified value of ${\ell }_r^{+}$. The zero is intentional: multiplying the minimum adjacent centroid step by k does not give a positive lower bound because the spatial route can turn. For distances between whole cells, let ${\rho }_r^{+}$ bound the distance from the reconstructed centroid to every point of any level r cell. Section 3.5.1 coverage derives such a bound from the base mesh. If distance between sets means their minimum point to point distance, then the distance from any ring cell footprint to the union of seed cell footprints lies in

$$\left[\max (0,L_k-2{\rho }_r^{+}),U_k\right].$$

(7)

The scale difference across levels is substantial. Near Calgary, the sampled adjacent centroid step at level 20 is approximately $2.5$ m. Four times that length is about 10 m. At level 28 the corresponding total is approximately $3.9$ cm. These are accumulated step scales, not minimum outward displacements: successive moves need not follow one geodesic. Neither these products nor $K=\lceil d/{\ell }_r^{+}\rceil$ guarantees coverage of a requested distance d. The inequality in Equation (6) gives an upper bound on reached distances, whereas completeness requires ruling out unreached cells within d.

The centroid distances in this subsection are evaluated after radial reconstruction on the sphere, not in a projected face plane. The gnomonic projection maps great circles to straight lines but is not distance preserving and has no single constant scale away from its center [36]; it therefore cannot by itself convert K into meters. Spherical offsetting similarly treats distance intrinsically on the sphere rather than as an uncorrected planar displacement [20]. All seed cells in one query must share a level. Mixed resolution neighborhood semantics and volumetric cells are outside the current scope.

3.3. Efficient Exact Halo Construction

The spherical interpretation shows why a ring count alone is not a distance guarantee. We first establish the reusable operation needed by both fixed hop and spatial distance queries: constructing complete rings from indexed seed cells. A cell representation must support persistent storage and hierarchy operations at the query interface, while allowing repeated adjacency and duplicate discovery with little work. This subsection introduces the two index representations and the local rules that connect them to the mesh.

Figure 5 shows the indexed halo component, not the complete spatial buffer procedure. Public hierarchical indices are converted once into compact integer triangular indices. Integer arithmetic constructs adjacent cells within a base mesh face, and an ACM record resolves a face boundary crossing. The CPU and GPU methods in Section 3.4 apply these rules to complete rings. Output indices are converted back after the requested expansion. The source cover and distance tests in Section 3.5 surround these calls and have their own cost.

This separation follows how often information is needed. Geographic geometry is used when creating the source cover and evaluating distances. Hierarchical decoding and encoding occur at the boundaries of a halo call. The repeated adjacency loop uses only integer coordinates, base mesh transitions, and the discovered set. Consequently, accelerating that loop does not remove the separate geometric work needed by a buffer query.

Figure 5. Architecture overview. Public hierarchical indices are decoded once into compact integer triangular indices. The repeated traversal reconstructs adjacent cell indices, consults an Atlas of Connectivity Maps (ACM) record only at a base mesh edge, and constructs exact rings through the discovered set. Returned cells are encoded once after all requested rings are complete.

Figure 5. Architecture overview. Public hierarchical indices are decoded once into compact integer triangular indices. The repeated traversal reconstructs adjacent cell indices, consults an Atlas of Connectivity Maps (ACM) record only at a base mesh edge, and constructs exact rings through the discovered set. Returned cells are encoded once after all requested rings are complete.

3.3.1. Public Hierarchical Index

The halo component receives cell indices rather than geographic coordinates. Its first requirement is a representation that identifies the input and output cells at the query interface. This representation must be persistent and reveal the base mesh face and resolution while supporting hierarchy operations before or after a halo query. The public hierarchical index is the unsigned 64 bit value accepted and returned by the query interface. Bit 63 is reserved and zero. Bits 62 through 57 store the base mesh face code. Bits 56 through 0 contain two bit child digits, followed by one level sentinel and trailing zeros. The prefix, sentinel, and padding convention, including recovery of the level from the least significant set bit, follows the mechanism documented for Google’s S2CellId [23]. Our adaptation assigns six bits to a triangular base mesh face code and supports 28 refinement digits. The base level sentinel is bit 56, and every refinement moves it two bits toward the least significant end. The 57 payload bits therefore hold 28 child digits and one sentinel.

The sentinel makes the level recoverable from the least significant set bit, and descendants of a cell share its prefix. Given a canonical level r public hierarchical index and a target level $t\le r$, an ancestor is computed by clearing the suffix below the target position and writing the target sentinel. The number of bit operations does not depend on $r-t$. The same prefix defines a descendant interval at a target level. We retain this representation at the query boundary because databases, files, ancestor queries, and workflows with several resolutions require the level and hierarchy to be recoverable from the index itself. Repeating this hierarchy processing for every same level cell adjacency would add instructions and memory traffic to the dominant query operation.

3.3.2. Compact Integer Triangular Index

The public hierarchical index preserves the refinement sequence, but the repeated adjacency operation needs different information: during a fixed level halo query, every frontier cell must generate three neighbors and participate in hashing with little repeated work. Carrying the complete refinement prefix through that repeated loop would preserve information that is not used until the result leaves the query. We therefore define the compact integer triangular index. It is a transient 63 bit index used only inside the fixed level expansion; $(f,i,j,o)$ are its coordinate fields, and i and j are the two integer coordinates within the chart of one base mesh face. Its layout is

$$\underset{6}{\underbrace{f}}\parallel \underset{28}{\underbrace{i}}\parallel \underset{28}{\underbrace{j}}\parallel \underset{1}{\underbrace{o}},$$

(8)

where $f\in [0,59]$ for the Pentakis case, i and j are integer chart coordinates, $o\in \{0,1\}$ is triangle orientation, and $N=2^r$. For upward cells $U(i,j)$, $i+j<N$; for downward cells $D(i,j)$, $i+j+1<N$. Figure 6 grounds these symbols in the subdivided chart of one base mesh face. Each 28 bit field stores integers from zero through $2^{28}-1$, so the layout supports $N=2^{28}$ subdivisions along an edge and sets the maximum supported resolution to $r=28$. The endpoint N is a geometric vertex coordinate, not a cell anchor stored in these fields. On a spherical Earth of radius $6,371.0088$ km, the $60·4^{28}$ Pentakis cells have a global average area of approximately $1.18$ cm². The square root of this average area is about $1.09$ cm, but it is not a triangle edge length or a certified cell diameter; individual areas and dimensions vary with projection distortion.

The relation to the continuous face coordinates is explicit. Let $({\beta }_0,{\beta }_1,{\beta }_2)$ be the barycentric coordinates of the projected point, with ${\beta }_0+{\beta }_1+{\beta }_2=1$, and set $N=2^r$, $x=N{\beta }_1$, and $y=N{\beta }_2$. For a point strictly inside one refined cell, its discrete index is obtained by

$$i=\lfloor x\rfloor ,j=\lfloor y\rfloor ,o=\left\{\begin{array}{cc}0,\hfill & \left\{x\right\}+\left\{y\right\}<1,\hfill \\ 1,\hfill & \left\{x\right\}+\left\{y\right\}>1,\hfill \end{array}\right.$$

(9)

Here $\left\{x\right\}=x-\lfloor x\rfloor$. Equation 9 is an interior rule, not a complete rule for boundary ties. At a point on a refined edge or vertex, assignment must instead select a valid incident cell whose closed footprint contains the point. One complete mathematical convention selects the smallest unsigned public hierarchical index among all such incident cells, including cells on neighboring base mesh faces. For example, $N=4$ and $(x,y)=(2,2)$ touch $U(1,2)$, $U(2,1)$, and $D(1,1)$ in the current face; the floor result $U(2,2)$ is invalid. At chart vertices A, B, and C, the valid incident cells in that face are $U(0,0)$, $U(N-1,0)$, and $U(0,N-1)$, respectively. This exact ownership convention is a mathematical specification. The measured floating point mapper uses face selection, endpoint clamping, and cell validation; the current evidence does not certify all exact boundary decisions against that convention. Conservative source covering is different: it retains every touching cell, as required in Section 3.5.1.

The ratios $i/N$ and $j/N$ therefore describe the grid vertex at integer coordinates $(i,j)$, whose barycentric coordinate is $(1-(i+j)/N,i/N,j/N)$; they are not, by themselves, the original continuous point coordinates. For completeness, the barycentric centroid of the indexed cell is

$$({\beta }_0^c,{\beta }_1^c,{\beta }_2^c)=\left\{\begin{array}{cc}\left(1-{\displaystyle \frac{i+j+2/3}{N}},{\displaystyle \frac{i+1/3}{N}},{\displaystyle \frac{j+1/3}{N}}\right),\hfill & o=0,\hfill \\ \left(1-{\displaystyle \frac{i+j+4/3}{N}},{\displaystyle \frac{i+2/3}{N}},{\displaystyle \frac{j+2/3}{N}}\right),\hfill & o=1.\hfill \end{array}\right.$$

(10)

The two orientation cases in Equation (10) reconstruct the planar cell centroid before radial normalization. After this discretization, the same level adjacency phase operates strictly on $(f,i,j,o)$; the continuous barycentric values are no longer required.

Child numbering before table acceleration.

The conversion tables accelerate a fixed refinement rule, so that rule must first define what each public child digit means. Let a parent have integer anchor $(i,j)$ and orientation o, and let $(i^{\prime },j^{\prime },o^{\prime })$ denote its child at the next level. A base mesh face begins at $(0,0,0)$. Digit 3 selects the central child and changes orientation; digits 0, 1, and 2 select the three corner children according to Figure 7 and Table 1. Iterating these rules gives the same compact integer triangular index as geometric midpoint subdivision.

Each public hierarchical index is decoded once before the first frontier. A read only table with 512 entries processes four refinement digits per lookup. Its nine bit lookup key contains the current orientation bit and the next four two bit child digits. Each entry returns the increments of i and j and the orientation after those four refinements. A level 28 public hierarchical index therefore requires seven lookups instead of 28 separate child steps. A second table with 512 entries performs the reverse conversion four levels at a time. It uses the current orientation and the lowest four bits of both i and j to recover four child digits and the parent orientation. The storage of both tables is fixed and does not depend on the resolution. Output cells are converted to public hierarchical indices only after all halo rings are complete. This placement avoids decoding a public hierarchical index for every generated candidate, while preserving the persistent representation needed by storage and hierarchy operations.

3.3.3. Comparison with BIM and Design Rationale

The preceding subsection introduced the compact integer triangular index for the repeated adjacency operation. This subsection documents its storage and operation requirements in relation to the more general Barycentric Indexing Method (BIM), rather than predicting a GPU speedup from coordinate counts. The two indexes share a barycentric foundation but serve different purposes and are not competing public DGGS standards. The published BIM index is $I=(i_b,w,a,b,c)$. It combines a base cell index $i_b$ and an explicit resolution w with the symmetric three axis integer coordinate $(a,b,c)$. BIM supports geometric, hierarchical, and neighborhood operations across watertight triangular meshes, and its Atlas of Connectivity Maps extension preserves these operations across base mesh boundaries [8]. Our compact integer triangular index is a temporary representation for one halo query at a shared resolution. The public hierarchical index remains the externally visible cell representation.

BIM retains three symmetric axes even though normalized barycentric coordinates have two independent degrees of freedom. As described in the BIM paper, its integer coordinate image is sparse, so some integer triplets do not represent cells. Point indexing near exceptional vertex configurations can therefore require a validity check and a small correction to select a represented cell [8]. Adjacency within one base triangle remains constant time because one of a, b, or c is incremented or decremented according to orientation. When an adjacency crosses a base mesh edge, the crossing axis and destination orientation determine one of three coordinate permutations: $(a,b)$, $(b,c)$, or $(c,a)$ is exchanged. This symmetry is valuable in a general mesh index. A direct tuple of three 28 bit coordinates occupies 84 coordinate bits, but that is not a lower bound on a packed BIM index. At a shared level, the permitted coordinate sums and a triangle type bit allow one coordinate to be reconstructed from the other two. The BIM paper explicitly discusses two coordinates plus one distinguishing bit. A specialized packed BIM key could therefore also fit two 28 bit fields, a 6 bit face field, and one type bit in 63 bits. Our choice is that the stored fields directly serve the implemented adjacency and hash operations, not that BIM inherently requires a multiple word key.

The compact integer triangular index omits state that does not change during this narrower workload. The level r is common to every seed, frontier, discovered, and output cell, so it is stored once for the query. After discretization, $(i,j)$ provide the two independent chart directions, and o distinguishes upward and downward elementary triangles. Every tuple that satisfies the relevant triangular domain condition represents one cell. Adjacency inside a base mesh face therefore uses the fixed offsets in Equations (11) and (12), without introducing a third coordinate into the interior operation. Valid cell traversal should not be confused with point assignment: this does not imply that BIM corrects an invalid triplet after every neighbor step. A base mesh edge crossing still requires a boundary decision, but one atlas record maps a scalar edge position using either $t^{\prime }=t$ or $t^{\prime }=N-1-t$. Together with the base mesh face code f, these fields form the 63 bit layout in Equation (8). Figure 8 the coordinate conventions and boundary operations used during a query.

Table 2. Coordinate conventions and execution roles. The comparison follows the published BIM operations; it claims neither an inherent packed storage advantage nor measured superiority over an optimized BIM implementation.

Design concern	BIM index	Compact integer triangular index	Halo query consequence
Coordinate image	Symmetric $(a,b,c)$ image is sparse; exceptional point assignment may require checks	Every tuple satisfying the explicit $U/D$ domain conditions represents a cell	Both methods distinguish point assignment from traversal of valid cells
State at level 28	A packed form can retain two coordinates and a triangle type, with face and query level	$f_6\parallel i_{28}\parallel j_{28}\parallel o_1$ uses 63 bits; r is query wide	Both can support a single word key; packing alone does not prove a speed advantage
Interior adjacency	Increment or decrement one of three axes according to orientation	Three fixed $(i,j)$ offsets and an orientation change	Two coordinate fields and no hierarchy decoding in the repeated loop
Base mesh boundary crossing	Detect the boundary axis and select one of three axis exchanges	Apply one of the $3F$ atlas records to scalar edge position t	Different coordinate rules; both need boundary handling and implementation dependent instructions
Intended role	General geometry, hierarchy, and movement operations	Adjacency expansion and set membership at one shared level	Generality remains at the public boundary; compact state serves the repeated query operation

Table 2. Coordinate conventions and execution roles. The comparison follows the published BIM operations; it claims neither an inherent packed storage advantage nor measured superiority over an optimized BIM implementation.

Design concern	BIM index	Compact integer triangular index	Halo query consequence
Coordinate image	Symmetric $(a,b,c)$ image is sparse; exceptional point assignment may require checks	Every tuple satisfying the explicit $U/D$ domain conditions represents a cell	Both methods distinguish point assignment from traversal of valid cells
State at level 28	A packed form can retain two coordinates and a triangle type, with face and query level	$f_6\parallel i_{28}\parallel j_{28}\parallel o_1$ uses 63 bits; r is query wide	Both can support a single word key; packing alone does not prove a speed advantage
Interior adjacency	Increment or decrement one of three axes according to orientation	Three fixed $(i,j)$ offsets and an orientation change	Two coordinate fields and no hierarchy decoding in the repeated loop
Base mesh boundary crossing	Detect the boundary axis and select one of three axis exchanges	Apply one of the $3F$ atlas records to scalar edge position t	Different coordinate rules; both need boundary handling and implementation dependent instructions
Intended role	General geometry, hierarchy, and movement operations	Adjacency expansion and set membership at one shared level	Generality remains at the public boundary; compact state serves the repeated query operation

Table 2 summarizes the design tradeoff. The complete compact integer triangular index fits one 64 bit discovery word, but a storage collision can still require several probes. BIM could use precomputed transformations or a different packing scheme; instruction counts and thread divergence depend on the implementation and compiler. We therefore make no measured GPU performance claim against an optimized BIM implementation. The conversion ablation in Section 4.1.2 instead compares keeping our compact integer triangular index throughout a query with repeatedly converting our public hierarchical index.

3.3.4. Implicit Local Adjacency

With the index roles established, we can derive the operation performed for every frontier cell. Inside one base mesh face, crossing an edge changes an upward triangle into a downward triangle, or conversely. The adjacent triangle either uses the same integer position or differs by one in one coordinate. These three possibilities give the fixed adjacency equations:

$$\begin{array}{cc}\hfill \Gamma \left(U\right(i,j\left)\right)& =\left\{D\right(i,j),D(i,j-1),D(i-1,j\left)\right\},\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill \Gamma \left(D\right(i,j\left)\right)& =\left\{U\right(i,j),U(i+1,j),U(i,j+1\left)\right\}.\hfill \end{array}$$

(12)

If a candidate satisfies the triangular domain conditions, it remains on face f and is encoded directly. Therefore, an adjacency inside one base mesh face performs no table lookup and loads no stored neighbor record. This reconstruction is the central scalability choice: storing three neighbor indices for every refined cell would require storage proportional to the complete level, whereas the arithmetic cost remains constant at every resolution.

3.3.5. Base Mesh Edge Atlas

The interior rules are sufficient only while the adjacent cell remains on the same face. At a base mesh edge, the query must identify the neighboring face and express the adjacent cell in that face’s coordinates. The ACM stores this relation once for each directed base mesh edge. For source face f and one of its three edges e, the record

$$T(f,e)=(g,e^{\prime },\rho )$$

(13)

stores the destination base mesh face g, destination base mesh edge $e^{\prime }$, and an edge direction reversal bit $\rho$. If $t\in [0,N-1]$ is the discrete position of the source boundary cell along edge e, then

$$t^{\prime }=\left\{\begin{array}{cc}t,\hfill & \rho =0,\hfill \\ N-1-t,\hfill & \rho =1.\hfill \end{array}\right.$$

(14)

The destination edge $e^{\prime }$ and mapped position $t^{\prime }$ select one of three fixed boundary templates, which reconstructs the destination compact integer triangular index $(g,i^{\prime },j^{\prime },o^{\prime })$. Thus, t and $t^{\prime }$ are local edge positions, not cell indices. In concrete terms, the ACM record answers the otherwise ambiguous boundary question: “which base mesh face owns the neighbor across this edge, which of its three edges is shared, and in which direction are the discrete edge positions ordered?” For a triangular base mesh with F faces, the atlas has exactly $3F$ directed records and is independent of refinement level. Pentakis has $60\times 3=180$ records and occupies 720 bytes in the packed representation; an icosahedron has 60 records and an octahedron has 24. The table resides in device constant memory because it is small, does not change during a query, and can be shared by threads crossing the same base mesh edge. Figure 9 contrasts the arithmetic only interior case with the exceptional boundary case that consults this record.

Table 3 completes this reconstruction. Edge 0 is the side opposite chart vertex A, edge 1 is opposite B, and edge 2 is opposite C, using the axes in Figure 6. Every cell sharing an edge with a base mesh boundary has upward orientation. For example, $N=4$, destination edge $e^{\prime }=0$, and mapped position $t^{\prime }=2$ reconstruct the destination cell $(g,1,2,0)$. Packing that tuple, not the scalar $t^{\prime }$, gives the destination compact integer triangular index.

The transition records must agree with the mesh geometry, not merely with another implementation using the same records. We generate the table from a directed edge representation, which records each face’s side of every base mesh edge. Independent references then check both the base transitions and complete refined edge pairs without consulting the runtime ACM. Section 3.6.6 gives the construction and test coverage. This separates the mechanism used to cross a face boundary from the evidence that the crossing is correct.

3.4. From Local Adjacency to Complete Rings

Section 3.3 established how one cell is represented and how its three edge adjacent cell indices are reconstructed. A complete halo query must coordinate these constant cost operations over every cell in the current ring while preserving the recurrence in Equation (1). At the start of iteration k, the frontier is the completed ring $H_{k-1}$ and the discovered set satisfies

$$V={\displaystyle \bigcup _{q=0}^{k-1}}{H}_q.$$

(15)

A candidate belongs to $H_k$ exactly when its compact integer triangular index is absent from V. This section first establishes the serial and multithreaded CPU realizations of that invariant. It then defines the GPU resident query state, the exact hash address and collision procedure, concurrent first discovery, synchronization between rings, capacity checks, and resource bounds. These steps answer the collective set construction question that remains after adjacency has been reconstructed for one cell.

3.4.1. Serial and Multithreaded CPU Traversal

The CPU implementations establish the ring recurrence and provide performance baselines with different storage and synchronization costs. All accept and return public hierarchical indices, use the same compact integer triangular indices and adjacency rules, and construct the same ring sets. Agreement checks execution consistency; independent geometry tests are still needed to validate their shared atlas.

The serial container implementation preloads the compact seed indices into a node based hash set and stores $H_0=P$ in the frontier. For every frontier cell it constructs three adjacent indices and appends an index to the next ring only if insertion reports that the index was absent. This implements Equation (15). The serial flat table baseline replaces that container with an array of open addressed slots. It uses the mixing and linear probing rule defined in Section 3.4.3 and retains table capacity between calls. This separates node allocation and pointer access from the cost of the recurrence.

Parallel processing also needs to coordinate insertions. The partitioned implementation uses the OpenMP shared memory programming interface to generate candidates, groups them by a hash selected shard through count and scatter, and gives each shard to one worker. The shard count is the next power of two above $\max (64,8T)$ for T threads, giving 256 shards at twenty threads. Each worker inserts into its shard without locks. New indices are concatenated after each ring. This design measures the cost of grouping and node based sets; it is not the primary multithreaded baseline.

The multithreaded flat table is the primary parallel CPU baseline. The primary matrix runs it with twenty threads, so the results refer to it as the twenty thread flat table; the thread sweep in Section 4.1.1 varies that count. It uses atomic compare and exchange to claim slots in a shared open addressed table. Each worker collects new cells in a private buffer; separately aligned headers reduce sharing of frequently updated counters. Dynamic scheduling distributes frontier work across the workstation’s performance and efficiency cores. Expansion and concatenation use separate OpenMP regions, with implicit synchronization at region boundaries; output encoding is also parallel. The CPU and GPU therefore share the membership and complete ring principles, but not an identical instruction sequence or barrier count. CPU output is materialized in a newly allocated returned container.

The flat table designs remove the node based set as the sole CPU reference. Section 4.1.1 compares both flat table baselines with the GPU and reports a thread count sweep. The serial container and partitioned implementations remain diagnostics of the baseline choice. The order of indices within a ring is unspecified; ring membership and ring offsets determine correctness.

3.4.2. GPU Resident State and Query Stages

A halo query is iterative: ring $H_k$ must be completed before it can serve as the frontier for constructing $H_{k+1}$. If each completed ring is returned to the CPU, every iteration requires GPU synchronization, transfer of intermediate state, a CPU decision about the next frontier, and another GPU launch. Because generating the three adjacent cells requires only a small amount of arithmetic, this repeated coordination can cost more than the adjacency operation itself. The purpose of the GPU resident design is therefore to retain the complete iterative state in GPU memory and execute all requested rings without transferring intermediate frontiers to the CPU. Figure 5 shows this placement: conversion occurs at the query boundaries, while adjacency expansion and discovery remain on the GPU.

Each GPU resident structure supports a specific part of the recurrence. The unique compact seed array initializes $H_0=P$. The discovered set V records the seed cells and every cell admitted to a halo ring, allowing later candidates to be rejected without consulting the CPU. The cumulative compact output stores the discovered halo cells and provides the frontier slices used by subsequent iterations. The $K+1$ ring offsets identify where each ring begins and ends, while device counters and an error status coordinate output allocation and report bounded failures. These structures remain in GPU memory throughout the expansion. Public hierarchical indices are used only in the input and output arrays and do not replace the compact integer triangular indices used by repeated adjacency and discovery operations.

The resulting query follows the four stages illustrated in Figure 5. First, the public hierarchical input indices are copied once to the GPU, validated, decoded, and deduplicated to form $H_0=P$ in compact form. Second, all requested ring iterations construct adjacent cell indices, test discovery, and append new cells while the discovered set, frontier slices, and ring offsets remain on the GPU. The final control state and completed ring offsets are then read by the CPU. Third, after the final ring is complete, the compact output indices are encoded into public hierarchical indices. Fourth, the encoded cell indices are copied once to the CPU. This separation also defines the reported timing components. Seed decoding and expansion share one kernel, but their recorded traversal interval also includes final control and ring offset readback, synchronization, and the intervening host scheduling interval before the end event is recorded. Output encoding and the bulk cell array transfer are measured separately; the traversal interval is not a pure kernel measurement. Geographic point projection is not part of this preindexed query.

3.4.3. Hash Function and Open Addressed Discovered Set

During a parallel frontier expansion, potentially thousands of GPU threads simultaneously construct adjacent cell indices. Different frontier cells can share the same adjacent cell, so several threads may generate an identical candidate. To preserve exact set semantics, the architecture must admit that cell to the new halo ring only once and reject it if it already belongs to the seed set or a completed ring. We support this decision using a GPU resident, open addressed hash table that stores the discovered set V. In an open addressed table, cell indices are stored directly in an array of slots. A lookup checks a sequence of slots until it finds the index or an empty location. Each slot check is called a probe. The following paragraphs explain how the sequence is chosen and how concurrent insertions are resolved.

Candidate construction produces one 63 bit compact integer triangular index x for each adjacent cell. The packed face and coordinate fields are spatially correlated, so directly using selected bits of x could place many nearby cells in the same part of the table. The mixing function $m\left(x\right)$ spreads these patterns before table addressing. It applies a fixed sequence of 64 bit addition, bitwise exclusive OR with shifted values, and multiplication by odd constants, following a standard SplitMix style integer mixing procedure [37]. These operations determine only the table location associated with x; they do not modify the compact integer triangular index itself.

Let the table capacity be $M=2^q$. Probe number ℓ addresses

$$p_{\ell }\left(x\right)=\left(m\left(x\right)+\ell \right)\&(M-1),0\le \ell <M,$$

(16)

where & denotes bitwise AND. Because M is a power of two, masking with $M-1$ maps the mixed value to a valid table slot. The first probe uses $p_0\left(x\right)$. If that slot contains a different index, linear probing advances through $p_1\left(x\right),p_2\left(x\right),\dots$ until the operation finds an empty slot, finds x, or exhausts the table [38].

Every valid compact integer triangular index has bit 63 equal to zero. The all ones word

$$E=2^{64}-1$$

can therefore represent an empty slot without conflicting with a valid cell index. The logical claim attempts to replace E with x atomically. The implementation first reads the slot and avoids the atomic operation when it is already occupied; Section 3.4.6 explains this optimization. When a claim is attempted, its returned value produces three outcomes:

1.

If the returned value is E, this thread inserted x and discovered a new halo cell.

2.

If the returned value is x, the index is already present, so the thread stops without appending it again.

3.

If the returned value is another index $y\ne x$, a storage collision occurred, so the thread advances to the next probe position.

The table is cleared before each query and preloaded with the unique seed indices. Every newly discovered halo index remains in the table for the rest of the query. No deletion occurs, so the table needs no additional marker for a deleted entry. Consequently, an index already present in the table may represent a seed cell, a cell discovered in an earlier ring, or the same candidate inserted by another thread during the current ring. All three cases require the same action: reject the candidate because it must not be appended again.

The target load

$$\alpha ={\displaystyle \frac{\left|V\right|}{M}}$$

is at most $0.5$. Under the standard assumption that the mixing function distributes indices across the table, this bounded load keeps the expected number of probes constant. Fewer probes reduce global memory accesses and contention among threads attempting to claim nearby slots. The explicit limit of M probes guarantees termination; exhausting all slots reports insufficient hash capacity instead of returning an incomplete halo. A different resident index y is therefore a hash storage collision, not a statement about spatial adjacency or halo membership.

Figure 10 connects this storage procedure to the spatial operation. It follows one adjacent cell from triangular candidate construction, through the 64 bit mixing function, to the atomic decision that appends the cell, rejects it as already present, or continues probing after a storage collision.

Figure 10. From triangular adjacency to atomic candidate discovery. A one into four refinement exposes the three edge adjacent candidates of a frontier cell $c\in H_{k-1}$. One candidate is represented by the compact integer triangular index x and passed through the 64 bit mixing function to select an initial probe slot. The atomic compare and swap returns Empty when x is first discovered, x when it is already present, or a different resident index y when a storage collision requires probing the next slot. Only the first outcome appends x to the new halo ring $H_k$.

Figure 10. From triangular adjacency to atomic candidate discovery. A one into four refinement exposes the three edge adjacent candidates of a frontier cell $c\in H_{k-1}$. One candidate is represented by the compact integer triangular index x and passed through the 64 bit mixing function to select an initial probe slot. The atomic compare and swap returns Empty when x is first discovered, x when it is already present, or a different resident index y when a storage collision requires probing the next slot. Only the first outcome appends x to the new halo ring $H_k$.

One table slot stores the complete compact integer triangular index, and one successful native 64 bit atomic compare and swap claims an empty slot for that index. A discovery may still require several probes and atomic attempts before it terminates. Output reservation is a separate atomic operation, described next. This distinction is important when interpreting the representation choice in Section 3.3.3.

3.4.4. Atomic First Discovery and Next Ring Construction

The hash procedure defines storage outcomes; the halo recurrence requires a set membership decision. This subsection establishes the correspondence between them. A lookup has two terminal membership results: the candidate is newly inserted or its compact integer triangular index is already present. The already present result has three possible query histories: the index belongs to the seed set, it was inserted in an earlier ring, or another worker inserted it during the current ring. Together with a successful insertion of a new cell, these are the four query meanings. The GPU resident hash table intentionally stores only membership, not the history that caused an existing index to be rejected.

The compare and swap that changes an empty slot from E to x is the instant at which the cell becomes discovered. This is the operation’s linearization point: although several workers may attempt discovery concurrently, only one can make that successful change. That worker alone may append x to the output for $H_k$. Workers that subsequently find x already present do not append it again.

Appending each accepted cell with a separate global counter update would create another shared point of contention. The implementation therefore combines appends within a warp, a group of 32 GPU threads, following the warp aggregated atomic technique documented by NVIDIA [39]. A ballot operation produces a bit mask identifying which threads successfully discovered cells. One thread reserves a consecutive output range for that group, and each successful thread uses its position among the set bits to select a unique location in the range. This reduces counter updates while preserving one append per discovered cell. Together, atomic discovery and grouped append construct the next ring without a separate candidate array, global sorting, and subsequent set subtraction.

3.4.5. Cooperative Multiple Ring Execution

Atomic insertion resolves individual discoveries, but a later ring must not begin until every cell in the current ring has been processed. This dependency requires synchronization across the workers, not a new distance or membership rule. GPU threads are organized into blocks, and the collection of blocks launched for a query is called a grid. The cooperative implementation limits the grid so that its blocks can remain resident together. Threads process successive frontier cells in strides, allowing that bounded grid to handle larger rings. A grid wide barrier waits for discovery to finish, one designated worker records the ring’s output boundary, and a second barrier makes that boundary visible before the next iteration. The new output slice then becomes the next frontier without returning its size to the CPU.

Algorithm 2 summarizes the device resident control flow. The notation Insert includes the mixing, probing, and atomic claim defined in Section 3.4.3; WarpAppend denotes the aggregated output reservation described above.

Algorithm 2 GPU resident halo traversal

Require:   public hierarchical indices for the input; shared level r; hop depth K 1: validate, decode, and deduplicate the input to obtain compact $H_0=P$ 2: preload P into discovered set V; initialize cumulative output count to zero 3: for $k\leftarrow 1$ to K do 4:     for all $x\in H_{k-1}$ in parallel do 5:         for all $y\in \Gamma \left(x\right)$ except the recorded discovering cell, when present do 6:            if $Insert(V,y)=New$ then 7:                $WarpAppend\left(y\right)$ to $H_k$ and record the edge from y back to x 8:            end if 9:         end for 10:     end for 11:     synchronize all workers 12:     one worker records the cumulative count as the end offset of $H_k$ 13:     synchronize all workers; expose the $H_k$ slice as the next frontier 14: end for 15: encode the complete compact output into public hierarchical indices 16: return public hierarchical indices and the $K+1$ ring offsets

The reported RTX 4070 measurements use this cooperative schedule. The implementation also synchronizes after one worker prepares the control state at the start of each iteration. Including the two barriers shown after expansion, it executes three grid wide synchronization points per hop and initialization synchronization before the loop. The primary fused decode and expansion kernel has two initialization barriers; the separately timed expansion kernel in the buffer study has one. Algorithm 2 expresses the ring dependency; this additional barrier makes its loop control safe across the grid. An automatic scheduling mode selects a host driven kernel per ring fallback when cooperative launch is unavailable; the fallback preserves the same ring sets but transfers control to the host between iterations. If the caller explicitly requests the cooperative schedule on unsupported hardware, the query reports that condition rather than silently changing the requested mode. When $K=0$, no ring expansion is performed and the returned halo is empty. In every schedule, the order of indices within one ring is unspecified.

3.4.6. Probe Reduction and Output Delivery

Keeping rings on the GPU removes intermediate host transfers, but not redundant discovery work or the cost of delivering a large result. Four mechanisms reduce those costs without changing the set recurrence. Their individual and combined effects are evaluated in Section 4.1.2.

Arrival edge skipping.

Every new cell was reached across a shared edge from a cell already discovered. The successful worker stores that reverse edge number in one byte beside the new output index. For an interior move, the orientation specific equations determine the reverse edge; for a base mesh crossing, the atlas gives the destination edge (). When the new cell is later expanded, this one neighbor is known to be discovered and need not be submitted for hash lookup again. The implementation still constructs all three geometric neighbors before skipping that lookup. Seeds have no recorded parent and submit three neighbors for lookup; later frontier cells submit two. The successful claim selects one parent, and ring synchronization completes its record before the next expansion.

Read before compare and swap.

An unsuccessful discovery often needs only to establish that an index is already present. The implementation uses a 64 bit global load that bypasses the first level cache to inspect a slot, and attempts atomic compare and swap only after reading the empty marker. An occupied slot never changes during a query. At the algorithmic level, an older empty value leads to an atomic attempt that resolves the current slot, while an occupied value permits direct comparison. This monotone state argument motivates the optimization. The current load is a cache controlled load of the Compute Unified Device Architecture (CUDA), not a C++ atomic load; output agreement tests are implementation evidence rather than a formal proof of device memory ordering.

Page locked output staging.

The query object retains a page locked (pinned) host buffer across calls. These pages are not paged out by the operating system during the transfer. Encoded output is copied directly into this buffer without allocating and initializing a new pageable result vector inside the readiness timer. The instrumented interval ends when the public indices are available in the pinned host buffer and the ring offsets have reached their separate host result vector. The query interface then copies the indices into its returned cell vector before returning to the caller. We therefore distinguish pinned host readiness from returned container completion; the former is not the complete call latency of the query interface.

Fused control initialization.

Counters and status words are grouped into one control allocation and cleared together. Free device memory is queried when buffers need to grow rather than on every call with sufficient retained capacity. Growth and capacity failure remain checked. This fourth mechanism bundles two reductions in fixed per call work; its ablation does not separate their individual effects.

These choices add one arrival byte per output capacity slot and retain host staging memory between calls. They are implementation techniques rather than new DGGS indexing principles. The timing boundaries in Table 6 distinguish their computation and delivery costs.

In the remainder of the paper, the GPU implementation denotes the cooperative implementation with all four mechanisms enabled. The four mechanism removal variant disables these options in the current source for a controlled comparison. We distinguish it from the earlier GPU release, which predates these mechanisms and supplies the historical schedule and global validation evidence. The E3 configuration and the separate WGS84 implementation also retain their own recorded build and timing scope (). Similar option settings do not establish identical source versions.

3.4.7. Capacity Planning and Explicit Failure

Keeping several rings inside one GPU launch also constrains memory allocation. All workers must continue using the same hash table and output addresses while the grid synchronizes, so these arrays are allocated before the launch rather than moved during it. The query therefore needs an explicit capacity estimate and a defined response when that estimate is insufficient. For the Pentakis implementation, the automatic output budget is

$$B_{\mathrm{auto}}=n_{\mathrm{in}}{\displaystyle \frac{3K(K+1)}{2}},$$

(17)

where $n_{\mathrm{in}}$ is the input count before deduplication. This expression follows the ring growth of the regular triangular cell adjacency and is used as the allocation budget for the Pentakis implementation. It is not asserted as a universal bound for every triangular base mesh. A different base mesh may supply a justified budget or use the topology independent degree three tree bound $3n_{\mathrm{in}}(2^K-1)$. The caller may also supply an explicit maximum output count. In every case, the capacity checks below protect result completeness if the supplied budget is insufficient.

For an output capacity B, the table capacity is

$$M=nextPow2\left(\max (1024,2(n_{\mathrm{in}}+B))\right),$$

(18)

which keeps occupancy at or below 50% when the result stays within B. Before buffers grow, a memory preflight rejects a request whose required allocation exceeds 85% of currently available device memory. During execution, separate status values report invalid or noncanonical input indices, mixed input levels, exhausted hash capacity, and exhausted output capacity. No partial halo is returned after any of these conditions. The caller may resubmit with a larger explicit bound or route the query to the CPU. Explicit failure and retry are used instead of during launch reallocation because moving the table or output array would invalidate addresses shared across the cooperative grid.

The recorded workspace metric can be reconstructed from these capacities. With M hash slots, two output arrays and an arrival byte array of capacity B, two input arrays of capacity $n_{\mathrm{in}}$, and $K+1$ ring offsets, its reported value is

$$A_{\mathrm{reported}}=8M+17B+16n_{\mathrm{in}}+8(K+1)+96\mathrm{bytes}.$$

(19)

This expression reproduces the 81 workspace values recorded in the primary matrix. The $17B$ term counts two 64 bit output arrays and one arrival byte per slot. It remains a reporting metric: the source allocates a 128 byte control block but counts 96 bytes, underreporting explicit device payload by 32 bytes at fresh capacities. Runtime state, static device data, allocator overhead, and the pinned host buffer are excluded. Retained allocations can also exceed the capacities requested by a subsequent smaller query. We preserve the recorded measurements and state their accounting boundary.

The automatic budget grows quadratically with depth even when overlap makes the realized halo much smaller. Its admissible input count therefore depends on the currently free memory and power of two table sizes, not a fixed fraction of nominal GPU capacity. The recorded boundary buffer experiment instead supplies an explicit output capacity of 200,000 cells to each call. Input storage and discovered set capacity remain separate; this is not a 200,000 cell limit on their combined occupancy. Exceeding the output capacity still causes explicit failure.

3.4.8. Complexity and Storage

We can now distinguish the work performed by the query from the memory reserved to execute it. This distinction matters because a conservative capacity estimate can exceed the number of cells actually discovered. Let $H={\bigcup }_{k=1}^K{H}_k$ denote the returned halo. Every expanded cell has three edge adjacent cells. To compare schedules consistently, we count the logical adjacency work before redundant parent edges are skipped:

$$W=3{\displaystyle \sum _{h=0}^{K-1}}\left|H_h\right|.$$

(20)

Under the bounded load and distribution assumption stated in Section 3.4.3, expected adjacency and discovery work is $\Theta \left(W\right)$. The current implementation still constructs those W geometric neighbors, but arrival edge skipping removes some before discovery; W is therefore not the number of lookup candidates, slot probes, or atomic instructions. Hash probing is not claimed to be constant in the deterministic worst case: one candidate may inspect at most M slots before reporting capacity exhaustion. Input decoding processes $n_{\mathrm{in}}$ supplied indices and costs $O\left(n_{\mathrm{in}}\right)$, even when duplicates reduce the unique seed set substantially. Clearing the hash table costs $\Theta \left(M\right)$, and output encoding and materialization cost $O\left(\right|H\left|\right)$ at the supported maximum level. The cooperative traversal also performs three grid wide synchronization points per requested ring, plus initial launch synchronization. These initialization and synchronization costs help explain why adjacency throughput alone does not predict the latency of a small query.

The base mesh topology occupies exactly $3F$ ACM records, or $\Theta \left(F\right)$ static storage. Occupied query state contains $\left|P\right|+\left|H\right|$ distinct cells, but allocated storage must also include unused capacity reserved before the cooperative launch. With M discovered set slots, $n_{\mathrm{in}}$ input indices, an output capacity B, and $K+1$ ring offsets, the actual allocation is

$$O(M+n_{\mathrm{in}}+B+K).$$

(21)

When the reserved capacities are within a constant factor of the realized unique input and output sizes, this reduces to $O\left(\right|P|+|H|+K)$. The automatic budget can reserve more memory when seed neighborhoods overlap or the input contains many duplicates; output sensitive occupancy does not by itself imply an equally tight allocation. In either case, storage is determined by the query and its chosen capacity, not by materializing the $F{4}^r$ cells or their refined adjacency records.

3.5. From Halo Rings to Spatial Buffers

The preceding method constructs exact rings from indexed seed cells, but it does not determine how many rings a distance request requires. A complete spatial buffer additionally requires a cover of the continuous source, a stopping condition that rules out undiscovered target cells, and a distance test applied to discovered cells. We now define these requirements independently of the execution method, first on the sphere and then in the numerical WGS84 boundary implementation.

3.5.1. Conservative Distance Coverage with Bounded Surplus

A request expressed in kilometers asks for complete spatial coverage, not a particular number of rings. In the Nose Hill Park example, every location in the 1 km buffer must be covered even where its boundary cuts through a cell. We now state the conditions needed to obtain that guarantee on the spherical DGGS, rather than only in the projected illustration. The basic idea is to continue exact ring expansion until an entire frontier is certified to lie beyond the requested distance, then retain all cells that could intersect the requested buffer. Uncertain boundary cells remain in the answer. Their maximum surplus distance is bounded by the cell extent, so precision can be chosen through the resolution without accepting omissions. This formulation follows the conservative inclusion principles of spherical offsetting and continuous coverage [20,21]; it adds a stopping certificate around our exact halo recurrence.

The region that must be covered.

Let $\mathcal{S}$ be a nonempty closed source geometry on the adopted sphere, and let $d\ge 0$ be the requested distance in meters. For the running example, $\mathcal{S}$ is the complete park boundary curve and $d=1000$ m. The definition also permits a filled polygon if its complete cell cover is supplied. Define

$$B_d\left(\mathcal{S}\right)=\{z\in {\mathbb{S}}_R^2:{dist}_{\mathbb{S}}(z,\mathcal{S})\le d\},T_{r,d}=\{c\in C_r:{\Omega }_c\cap B_d\left(\mathcal{S}\right)\ne ⌀\},$$

(22)

where ${\Omega }_c$ is the closed spherical footprint of cell c and distance to $\mathcal{S}$ means the shortest great circle distance to any point of that geometry. $T_{r,d}$ contains all and only the cells that intersect the buffer, including its boundary. We call this the exact cell intersection set. A conservative returned set $Q_{r,d}$ must contain $T_{r,d}$, but may include cells that do not intersect the buffer. We also require its complete footprints to remain inside $B_{d+\eta }\left(\mathcal{S}\right)$ for a stated surplus width $\eta$. These are separate requirements: no missing cells and a bounded outer extent.

The seed set must include every cell whose closed footprint intersects $\mathcal{S}$. In particular, a point or curve on a cell edge or vertex contributes all touching cells; the unique ownership rule used for ordinary point indexing is not sufficient for this cover. The same formulation can accept other source geometries. A filled polygon includes its interior. A boundary curve, such as a park boundary or a shoreline, instead produces a buffer on both sides of the curve, and all components and islands must be included. A filled polygon must not silently replace a boundary curve.

A conservative cell distance test.

We need a test that cannot discard a cell merely because its centroid lies just beyond d. Let $s_c$ be its reconstructed spherical centroid and let ${\rho }_c$ be a certified upper bound for every distance from $s_c$ to a point of ${\Omega }_c$. Suppose a geometric routine supplies an interval $[a_c,b_c]$ containing ${dist}_{\mathbb{S}}(s_c,\mathcal{S})$, with width at most $\epsilon$. The triangle inequality gives a lower bound for distance from the whole cell to the source:

$$\lambda \left(c\right)=\max (0,a_c-{\rho }_c)\le {dist}_{\mathbb{S}}({\Omega }_c,\mathcal{S}).$$

(23)

The conservative acceptance test is $\lambda \left(c\right)\le d$. Equality and unresolved numerical cases are retained. A comparison of the centroid distance alone with d would miss cells that straddle the requested boundary. The interval routine must bound distance to the complete source geometry, not only its sampled points; outward rounding must also be used when forming the test and its radius bound.

Certifying how many rings are sufficient.

We keep the topological recurrence unchanged and write $V_K={\bigcup }_{q=0}^K{H}_q$ for all discovered cells, including the seeds. Starting from the complete source cover, expand rings until the first tested $K\ge 1$ for which

$$H_K=⌀\mathrm{or}\underset{c\in H_K}{\min }\lambda \left(c\right)>d.$$

(24)

The resulting depth is denoted by $K_{★}$; the minimum is evaluated only for a nonempty frontier. It is taken over the entire ring, across every source component and every base mesh face. A camera view, a few boundary cells, or a mean step length cannot replace this test. Intermediate rings are not pruned by the distance test, so their original minimum topological distances remain unchanged. If a valid upper estimate of K is available it can schedule work, but acceptance of the result still requires the certificate.

Proposition 1  

(Sufficiency of the whole ring stopping test). Assume a finite, conforming triangular tessellation of the sphere with closed cell footprints and connected incident cell fans. Let $\mathcal{S}$ be nonempty and closed, let P contain every cell touching $\mathcal{S}$, and let every $\lambda \left(c\right)$ be a valid lower bound on distance from the whole footprint to $\mathcal{S}$. If the complete rings satisfy Equation (24) at depth K, then every cell in $T_{r,d}$ belongs to $V_K$.

Proof. 

For a cell in $T_{r,d}$, select a point where its footprint meets $B_d\left(\mathcal{S}\right)$. A shortest spherical arc to a nearest point of $\mathcal{S}$ exists and remains in $B_d\left(\mathcal{S}\right)$. Cells intersected by this arc connect to a seed through edge adjacency. At an edge or vertex contact, include all touching closed cells; the connected incident fan supplies an edge adjacent route, and every included cell still intersects the buffer. If the target were outside $V_K$, this route from a seed would cross $H_K$, because minimum adjacency depth changes by at most one across an edge. That crossing cell has $\lambda \left(c\right)\le d$, contradicting the nonempty ring condition in Equation (24). If $H_K$ is empty, no further cell is reachable from the seeds; the same route again excludes an undiscovered target.    □

Proposition 1 is a geometric statement under explicit assumptions, not a proof that floating point arithmetic produces the required bounds. The WGS84 extension must satisfy the corresponding assumptions in its own metric; its numerical checks are reported separately.

After certification, retain

$$Q_{r,d}=\{c\in V_{K_{★}}:\lambda \left(c\right)\le d\}.$$

(25)

The seed ring $H_0$ is included when selecting this distance result, although the topological halo H excludes it. Otherwise, a requested location inside a seed cell could be omitted. The complete procedure is stated in 3.

Algorithm 3 Conservative distance cover using complete halo rings

Require:  closed source $\mathcal{S}$; distance d; level r; certified distance and cell radius bounds 1: $P\leftarrow$ a cover containing every cell touching $\mathcal{S}$ 2: $H_0\leftarrow P$; $V\leftarrow P$; $k\leftarrow 0$ 3: repeat 4:     $k\leftarrow k+1$; $H_k\leftarrow \Gamma \left(H_{k-1}\right)\setminus V$ 5:     $V\leftarrow V\cup H_k$ 6:     evaluate $\lambda \left(c\right)$ for every $c\in H_k$ 7:     fail explicitly if capacity, arithmetic, or certification requirements cannot be met 8: until $H_k=⌀$ or ${\min }_{c\in H_k}\lambda \left(c\right)>d$ 9: evaluate missing $\lambda \left(c\right)$ values in V, including seeds 10: return $Q_{r,d}=\{c\in V:\lambda \left(c\right)\le d\}$, $K_{★}=k$, and surplus bound $\eta$

Bounding additional coverage.

The same cell radius that prevents omissions bounds the excess. For an accepted cell, $a_c\le d+{\rho }_c$ and its true centroid distance is at most $a_c+\epsilon$. Every point in that cell is therefore at most $d+2{\rho }_c+\epsilon$ from $\mathcal{S}$. Assuming complete source covering and valid distance intervals, with ${\rho }_r^{+}\ge {\max }_c{\rho }_c$, the mathematical result satisfies

$$T_{r,d}\subseteq Q_{r,d},B_d\left(\mathcal{S}\right)\subseteq \bigcup _{c\in Q_{r,d}}{\Omega }_c\subseteq B_{d+\eta }\left(\mathcal{S}\right),\eta =2{\rho }_r^{+}+\epsilon .$$

(26)

This bounds the additional distance, not the percentage of extra cells or area. Those percentages depend on the length and shape of the source and on the land distribution around it. A bound on $K{\ell }_r^{+}$ alone provides neither inclusion in Equation (26).

The cell radius can be bounded without enumerating the refined grid. For base mesh face f, let $g_f$ be its planar centroid, $h_f>0$ its plane distance from the sphere center, and $m_f$ the maximum Euclidean distance from $g_f$ to its three vertices. Midpoint subdivision scales every planar cell by $2^{-r}$. Along a segment in that plane, radial normalization changes angle by at most $1/h_f$ per unit planar distance. Any radius value satisfying

$${\rho }_r^{+}\ge {\displaystyle \frac{R}{2^r}}\underset{f}{\max }{\displaystyle \frac{m_f}{h_f}}$$

(27)

bounds the distance from each normalized planar centroid to every point of its spherical cell. The maximum involves only base mesh faces. Adjacent cells share a point, so the same argument gives ${\ell }_r^{+}\le 2{\rho }_r^{+}$ and therefore $U_k\le 2k{\rho }_r^{+}$, still only an outward reach bound. The accompanying exact rational check verifies $m_f/h_f<0.44$ for all 60 faces of the frozen Pentakis mesh, a serialized and checksummed copy of the base mesh. Using ${\rho }_r^{+}=0.44R/2^r$ with $R=6,371,008.8$ m gives the rounded upward bounds in Table 4. This check certifies the mesh constant relative to its serialized coordinates; it does not certify an untested geographic distance routine or establish that every performance experiment used a byte identical mesh file.

The guarantee is relative to the supplied geometry, the chosen Earth model, and certified arithmetic. A spherical 10 km distance is not identical to a World Geodetic System 1984 ellipsoidal distance. An ellipsoidal contract requires compatible geodesic routines and corresponding bounds [35]. An incomplete source dataset or an unbounded numerical approximation cannot support a claim of complete real world coverage.

The preindexed traversal measurements evaluate the fixed K recurrence. The GPU returns complete topological rings and does not itself trim them against the spatial buffer. The separate extension below adds source covering, whole ring stopping tests, and cell filtering. Unknown $K_{★}$ requires a capacity limit and either resumable expansion or safe replay; a capacity error does not permit an incomplete answer. Cell filtering examines all visited cells, including seeds, rather than only the final ring. Even the exact cell intersection set generally extends beyond a continuous distance contour when its whole cells are displayed. Matching that contour exactly would require clipped cell geometry.

3.5.2. WGS84 Boundary Buffers and Numerical Implementation

The spherical formulation explains why whole ring tests can establish coverage. A request measured on WGS84 needs the same reasoning with an ellipsoidal metric and the actual curved cell footprints. We therefore retain the same ring recurrence and add an ellipsoidal source cover, stopping test, and cell distance filter. Its source $\mathcal{S}$ is the union of the supplied shortest WGS84 geodesic segments. Polygon rings contribute their boundaries, not their filled interiors. All components are preserved. This matches the boundary request introduced in Figure 1, including any central region that is excluded because it lies farther than d from the curve. Interpreting the input segments as geodesics is explicit because generic GeoJSON edges are not defined that way.

The WGS84 model has equatorial radius $a=6,378,137$ m and flattening $f_E=1/298.257223563$. Direct and inverse geodesic calculations use GeographicLib and the algorithms of Karney [35]. The spherical DGGS longitude and latitude become geodetic longitude and latitude on WGS84. This mapping preserves the cell indices and adjacency. It does not replace cell edges by shortest ellipsoidal segments: each footprint is the image of the complete radially projected planar triangle.

To bound these footprints, let $b=a(1-f_E)$ and let $M\left(\varphi \right)$ and $N\left(\varphi \right)$ denote the meridional and prime vertical radii of curvature at geodetic latitude $\varphi$. Both are bounded above by $a^2/b<6,400,000$ m. The ellipsoidal line element is therefore at most $6,400,000$ times the corresponding unit sphere line element. For any planar triangular patch $\tau$ inside base mesh face f, let $m_{\tau }$ be its greatest vertex distance from its planar centroid. The same radial derivative argument as in Equation (27) gives

$${\rho }_{\tau ,E}=6,400,000{\displaystyle \frac{m_{\tau }}{h_f}}$$

(28)

as a geometric radius bound in meters before numerical allowances. This applies to the actual curved footprint, including interior points. The verified base ratio $0.44$ gives rounded upward whole cell bounds of $2.685547$ m at level 20 and $0.010491$ m at level 28. The spherical values in Table 4 must not be substituted for these ellipsoidal values.

Distance to one source segment.

The cell test needs the distance from a point to the complete source curve, not only its vertices. Let $\gamma :[0,s]\to \mathcal{S}$ be one segment parameterized by geodesic arc length. Let $D_0$ and $D_1$ be the WGS84 distances from a query point q to its endpoints. Before numerical allowances, the triangle inequality gives

$$a_{\gamma }=\max \left(0,{\displaystyle \frac{D_0+D_1-s}{2}},D_0-s,D_1-s\right),b_{\gamma }=\min (D_0,D_1).$$

(29)

The endpoint distances give an upper bound because those endpoints belong to the segment. The lower bound follows because distance to q can change by no more than the distance traveled along the segment. Evaluating its geodesic midpoint creates two shorter intervals with their own endpoint bounds. The routine repeatedly divides the interval with the smallest lower bound, retaining the best sampled upper bound, until their difference meets the requested tolerance.

For a synthetic arithmetic example, a segment of length 10 m with endpoint distances both 13 m gives $[8,13]$ m in Equation (29). If its midpoint distance is 12 m, the combined interval tightens to $[10,12]$ m. This explains interval refinement; it is not a measured DGGS query.

Searching a complete boundary efficiently.

Most segments of a detailed boundary cannot contain the nearest point. A segment of length s is enclosed by a metric ball centered at its geodesic midpoint with radius $s/2$. Groups of segments are enclosed by larger balls: for a chosen group center, the radius is the maximum distance to a member segment midpoint plus that segment’s radius. The triangle inequality guarantees containment. Distance to the group center minus its radius supplies a lower bound. A group is skipped only when this bound exceeds an already established upper bound on the nearest source distance. The hierarchy reduces search work without replacing the source by samples.

Covering the complete source.

Source covering starts with every base mesh face. The point to source interval and patch radius give a lower bound on the distance from each closed triangular patch to the source. A patch is discarded only when that lower bound is strictly positive; remaining patches are refined to the query resolution. Equality is retained. A source touching a shared edge or vertex is therefore covered by every incident touching cell, rather than assigned to one owner. This statement assumes valid numerical bounds. It explains why the point ownership convention used for ordinary indexing is not used for source covering.

Stopping expansion and filtering cells.

The host loop requests complete halo queries at depths $1,2,4,\dots$ until the entire final ring has footprint distance lower bounds greater than d. Each request replays the recurrence from the same seed set. Iterations within a GPU halo call remain on the device, but distance tests between calls run on the CPU. The complete boundary workflow is therefore hybrid.

Final filtering examines every discovered cell, including the seeds. A lower bound greater than d excludes a footprint; an upper bound at or below d confirms intersection. For an uncertain cell, subdivision tightens the bounds over its curved footprint. This subdivision is internal to the geometric test: it does not change the cell’s public hierarchical index or the query’s output resolution. Cells still uncertain when a stopping threshold or evaluation budget is reached are retained and labeled. A finer query resolution instead changes the returned cell size and is a separate choice.

Numerical settings and their scope.

The implementation uses ordinary double precision and three distinct numerical settings. A $10\mu$m engineering allowance enlarges distance intervals and footprint radii. For a patch of radius ${\rho }_{\tau }$, the segment routine uses target interval width ${\tau }_{\mathrm{seg}}=\max (0.001\mathrm{m},0.01{\rho }_{\tau })$. The cell classifier stops as unresolved when its footprint interval is no wider than $0.0001$ m or its budget of 4096 evaluated patches is exhausted. Because the segment tolerance is not reduced below 1 mm, the classifier threshold does not imply uniform 0.1 mm distance accuracy.

The actual interval and classification status determine whether a cell is confirmed, excluded, or retained as unresolved. In Equation (26), $\epsilon$ bounds the width of the distance interval actually supplied to the test; it is not automatically the classifier’s stopping threshold. These are inspectable engineering controls, not outward rounded interval arithmetic. Zero unresolved cells means classification completed under this numerical contract, not a formal floating point certificate or proof of survey accuracy. Insufficient source, output, distance evaluation, or hop capacity reports an error rather than an accepted partial buffer.

3.6. Experimental Design

The evaluation addresses four questions in sequence. First, is GPU halo construction faster than serial and multithreaded CPU execution when all implementations use the same sparse representation, and for which workloads? Second, which representation and execution choices contribute to the gain, and how does the workload affect it? Third, what changes when geographic mapping, source covering, and output materialization are charged to the query? Finally, do independent geometric references and execution tests support exact adjacency and ring membership? Correctness checks are prerequisites for accepting every timing result even though their detailed results are reported last. Complete WGS84 boundary queries use a separate experiment and timing contract.

The following subsections define the real world datasets, global validation inputs, query configurations, hardware platform, CPU and GPU baselines, timing boundaries, component ablations, and acceptance criteria used to answer these questions.

3.6.1. Real World Datasets

The overlap and arrangement of seed cells affect how many candidates are generated and how often the same cell is encountered. We therefore select municipal datasets with three different spatial structures: clustered buildings, connected street lines, and an extended boundary curve.

We use three official City of Calgary datasets: Buildings roof outline polygons [40], Street Centreline [41], and City Boundary [42]. Acquisition uses the official Socrata data service with stable row identifier pagination, remote feature counts, nonempty geometry checks, atomic file writes, and Secure Hash Algorithm 256 (SHA 256) checksums. The downloaded sources contain 511,586 buildings, 120,008 street features, and one municipal multipolygon.

The datasets are converted into deterministic point distributions in the projected coordinate reference system with European Petroleum Survey Group (EPSG) code 3776 for metric sampling and then returned to World Geodetic System 1984 longitude and latitude. Buildings contribute 100,000 representative points from selected footprints. Street lines are sampled at 20 m spacing in source order, capped at 250,000 raw samples, and reduced to 223,008 distinct points. The municipal boundary is sampled at 10 m spacing to 19,717 points. Fifteen invalid building geometries are repaired during preparation; no street or boundary geometry requires repair. Table 5 reports the source and prepared counts, while Figure 11 shows the three resulting spatial distributions. These inputs test traversal structure, not polygon cover accuracy.

3.6.2. Experiment Matrix

To separate source size, refinement level, and hop depth, we vary each factor independently and repeat the same cells across the serial CPU, multithreaded CPU, and GPU implementations. This factorial structure supports comparisons of complete preindexed halo calls and controlled component ablations; the boundary buffer experiment uses a separate matrix.

For each dataset we request evenly strided subsets of 1,000, 10,000, and 100,000 points from the sorted comma separated value (CSV) file. A request greater than the available boundary samples uses all 19,717 points. Point to cell mapping removes duplicate cells before traversal. We test levels 12, 20, and 28 and hop depths 1, 4, and 8. The full factorial design is

$$3\mathrm{datasets}\times 3\mathrm{sizes}\times 3\mathrm{levels}\times 3\mathrm{depths}=81$$

configurations per measurement matrix. Each configuration uses 10 warmup calls and 30 timed repetitions.

Hypotheses behind the controlled variants.

One variant changes the execution domain without changing the DGGS cells. A rhombic macrochart groups two adjacent base mesh faces into one compute domain, so their common base mesh edge can be handled internally. Pairing the 60 triangles produces 30 such domains. The hypothesis is that fewer between domain transitions improve locality; the test compares natural triangles and macrocharts under the same sorting procedure, so it does not attribute a sorting cost to geometry alone.

The controlled arrangements also isolate spatial overlap. Spatially compact seeds are a deterministic prefix of a growing topological disk, whereas scattered seeds are distinct cells generated using a fixed pseudorandom sequence of face and child digits. Both use the same compact integer triangular index. Here spatial compactness describes arrangement, not index packing. Boundary centered seeds stress crossings between base mesh faces.

The primary 81 query matrix compares four CPU implementations with a GPU implementation using NVIDIA’s CUDA platform. A separate level 28 component matrix uses six controlled workloads: 10,000 spatially compact cells at four hops; 100,000 spatially compact cells at one, four, and eight hops; 100,000 scattered cells at four hops; and 100,000 cells centered on base mesh boundaries at four hops. It compares (i) host managed sets and hop control, (ii) GPU expansion followed by sorting and duplicate removal after every hop on natural triangles, (iii) the same sorting method with rhombic macrocharts, which group adjacent base mesh faces into larger processing domains, and (iv) the earlier GPU release, before the mechanisms of Section 3.4.6 were added. This schedule study was run on the release before the four mechanisms were added and has its own reference and timing contract. A dual index study separately expands 100,000 spatially compact seeds through four hops at levels 16 and 28. It compares retaining the compact integer triangular index throughout the query with a public hierarchical index only variant that decodes and encodes the hierarchy during every expansion. Every variant must return the same materialized set. Separate ablation sessions change one mechanism or one architectural choice at a time on five representative queries, as described in Section 4.1.2.

The dual index study uses a reproducible compact input generated in memory rather than a saved seed file. Its center lies on base mesh face 17 and follows the repeating child sequence $1,0,3,2$ to the requested level. The generator constructs a disk through radius 317 and retains the first 100,000 cells, in ring order and then first discovery order. For $U(i,j)$, neighbors are visited as $D(i,j)$, $D(i,j-1)$, and $D(i-1,j)$; for $D(i,j)$, the order is $U(i,j)$, $U(i+1,j)$, and $U(i,j+1)$. The retained input has exactly 100,000 unique cells: complete rings through radius 257 and 540 cells of ring 258. At levels 16 and 28, respectively, the center is $(17,38550,6939,0)$ and $(17,157903209,28422577,1)$. Each four hop comparison uses five warmups and thirty timed repetitions. The two variants must return the same cell set and checksum at each level.

The levels deliberately cover coarse, fine, and ultra high resolution indexing. Using a spherical Earth area and the global cell count, their mean cell areas are approximately $5.07\times {10}^5$ m², $7.73$ m², and $1.18\times {10}^{-4}$ m² at levels 12, 20, and 28, respectively. Fixed hop counts are therefore computational workloads, not matched metric buffers.

Spatial scale and planar checks.

The distance study samples 2,048 indexed cells per dataset and measures each cell’s three adjacent centroid distances, giving 6,144 edge samples per dataset at each tested level. We report medians and percentile intervals, rather than interpreting a sampled step length as a coverage guarantee. Separate planar checks compare exact point to triangle intersection distances with the conservative cell radius test. A reference case uses 1.5 km triangle sides and a 10 km request. A further factorial check uses side lengths $0.8$, $1.5$, and $2.5$ km, requested distances 0, 1, 4, and 10 km, and source coordinates $(0,0)$ and $(0.27,0.19)$ km, giving 24 configurations. The two source positions lie at a grid vertex and inside a cell. We also enumerate intersections with the projected 1 km Nose Hill boundary buffer used in Figure 1. The point source checks use a ${10}^{-10}$ km acceptance tolerance to test omissions and the surplus envelope. They are double precision illustration checks, not certified interval arithmetic, spherical validation, or GPU timing experiments.

3.6.3. System and Implementations

Launch cost, memory access, and CPU coordination depend on the platform. Measurements use Windows 11, an Intel Core i7 12700K with eight performance cores and four efficiency cores (twenty logical processors), and an NVIDIA GeForce RTX 4070 with approximately 12 GB of device memory and compute capability 8.9. The build used for the primary matrix and the ablation sessions records CUDA toolkit 13.3.73 and Microsoft Visual C++ (MSVC) 19.51.36248.0 with C++17. The release used for the topology validation and the schedule study recorded CUDA 13.2.51, and the separate boundary buffer study used CUDA 13.3. These are distinct experiment versions, not interchangeable descriptions of one build. Driver 595.79 is documented in the manifests of that release; the timing harness of the primary matrix does not itself capture a fresh driver record.

The primary GPU configuration enables the four mechanisms of Section 3.4.6 and disables scatter clearing. These switches must be supplied explicitly because the source defaults do not enable them. A measurement session is a named collection of measurements with recorded configurations and settings. A primary matrix session covers all 81 configurations, each in a fresh process that executes the included implementations in a fixed CPU before GPU order. Four sessions contribute to the primary matrix and are identified by the labels of their raw files: E3 ran the four CPU implementations and its GPU configuration without the four mechanisms, ALLOPT ran the serial flat table, partitioned implementation, twenty thread flat table, and a GPU configuration with scatter clearing enabled in addition to the four mechanisms, and A4s1 and A4 ran the serial flat table, the twenty thread flat table, and the GPU implementation. Two further sessions, E2B (removal of one mechanism at a time) and E2C (one architectural alternative at a time), supply the controlled variants of Section 4.1.2; E1, E2, and SCALE are earlier stage, ablation, and thread sweep sessions retained in the released files. Primary matrix measurements use the same mapped cells, ten warmups, and thirty timed calls per included implementation. The CPU thread sweep separately tests one through twenty workers. No core affinity policy is established by the saved records, and selecting twenty workers does not by itself establish full utilization or which core type performs each operation.

3.6.4. Timing, Metrics, and Acceptance

Performance depends on where a timer starts and how the result is delivered. The primary GPU interval measures pinned host readiness: input transfer, discovered table initialization, seed decoding and ring expansion, output encoding, and transfer into the retained pinned host buffer. It stops before the result is copied into the returned vector of the query interface. CPU intervals include construction of their returned vectors. The primary ratio therefore compares different storage ownership boundaries; an additional GPU container copy estimate is reported beside the largest result. Neither readiness nor the sum of recorded stages is asserted to be a wall clock timer around the complete call.

Warmup reuses capacity and removes initial setup from steady state measurements. Point mapping is timed once per process after file reading; it includes point assignment, sorting, and duplicate removal. Downloading, repairing source geometry, preparing samples, and file input and output are excluded. Mapping inclusive values add this measured preprocessing cost to each implementation’s halo summary; they are reconstructed timing views, not separately repeated complete workflows.

Each measurement summary records latency statistics, counts, workspace, transfers, and a result hash. The measurement harness retains the last output of each implementation and compares its complete sorted cumulative set against the reference. This is not a separate verification of all thirty calls or of every ring. The global ring tests and independent geometry tests have the distinct scope described in Sections 3.6.6 and 4.3.

Session selection and sensitivity.

The released matrix draws CPU records from four full sessions and the GPU implementation records from two of those sessions. The GPU records of E3 and ALLOPT use their respective alternative GPU configurations and do not add repetitions of the GPU implementation. Each implementation summary is eligible when the maximum of its thirty calls is at most three times the minimum (the call spread criterion). At least one eligible summary exists for every implementation and configuration. The reported latency is the smallest eligible session median for that implementation and configuration. This is a best observed session summary, not evidence that the chosen session was disturbance free.

For each configuration the comparison baseline is the smaller selected median of the serial and twenty thread flat tables. The ratio is this CPU value divided by the selected GPU readiness median. Selection is independent across implementations: 73 of 81 ratios use CPU and GPU records from different sessions. Each process executes its included implementations in a fixed order, but independently selected ratios are not paired session observations. Taking minima also does not prove an unbiased ratio or that GPU gains cannot be inflated. We therefore report matched session sensitivity in Section 4.1.1. The fixed execution order and unequal numbers of eligible sessions remain limitations.

Stage decomposition uses directly recorded means from the mechanism removal session (E2B), not rescaled estimates of median stages. The mechanism removal and architectural variant sessions each have their own GPU implementation reference, and each variant runs in a separate process. Ratios use the reference from the same experiment stage. They are not simultaneous paired runs and have no measured confidence intervals. The combined four mechanism removal is explicitly distinguished from changing one option.

3.6.5. Routing Between Implementations and Internal Checks

Small queries may finish on one CPU thread before GPU launch and delivery costs are amortized. A simple local selector uses two quantities known before traversal: level r and the automatic output budget $B_{\mathrm{auto}}$ from Equation (17). The proposed rule chooses GPU execution when $r\ge 20$ and $B_{\mathrm{auto}}\ge 5,000$; otherwise it chooses the serial flat table. Memory feasibility is checked separately. This is an empirically evaluated rule for the sampled Calgary workloads, not a topological theorem or a demonstrated dispatcher for arbitrary input extents.

An internal check withholds all 27 configurations of one dataset and fits only the bound threshold on the other 54, keeping $r\ge 20$ fixed. The threshold grid is $\{1000,2000,3000,4000,5000,6000,8000,10000,15000,20000,25000,30000,40000,50000\}$. The smallest threshold attaining maximum training accuracy is selected. Predictions are then tested on the withheld dataset. The target compares GPU readiness with serial flat table latency; the selected records also permit a check against the faster of both CPU flat tables.

This split withholds a dataset, not a city, hardware profile, level, or depth. The rule family and level condition were motivated by the same matrix, so the test is an internal threshold check rather than independent validation of the complete design. Input is deduplicated before these experiments, making supplied count equal to $\left|P\right|$. Duplicate input or a broader geographic extent can invalidate the observed relationship between level, capacity bound, and realized output. A source extent cell count is not used as a proven upper bound: halos can grow beyond that extent.

3.6.6. Independent Global Topology Validation

The most important correctness question is whether the Atlas of Connectivity Maps represents the actual oriented base mesh. Equality between CPU and GPU outputs cannot answer this question because both implementations use the same atlas records. The retained validation evidence has three kinds; the first two are independent geometry references, while the third compares complete CPU and GPU ring sets. These tests were run on the earlier GPU release and have not been repeated with all four mechanisms enabled. First, a frozen oriented half edge representation of the base mesh, serialized independently of the runtime atlas, reconstructs all 180 expected directed base transitions without reading that atlas. Exhaustive adjacency evaluation over levels 0 through 5 checks 81,900 cells, 245,700 directed adjacency operations, and 11,340 directed positions on base mesh boundaries. Second, a level 6 reference constructs canonical integer indices for refined vertices and pairs complete refined edges without using atlas transitions. Its 245,760 cells contain 368,640 undirected edges, yielding 737,280 independently derived directed adjacency relations and 11,520 directed adjacencies across base mesh boundaries.

Third, a global query matrix distributes seeds over every directed base mesh edge and around every base mesh vertex. The edge scenario samples five positions on each of the 180 directed edges. The vertex scenario covers the 12 Pentakis vertices with valence five and the 20 vertices with valence six through their incident charts. We test levels $\{2,6,12,20,28\}$ and depths $\{1,2,4,8\}$, producing 40 CPU and GPU comparisons. Every ring is sorted and compared, and the matrix records 56,160 base mesh edge crossings. A separate performance case places 100,000 level 28 sources on base mesh boundaries, expands them through eight hops, and records 1,447 crossings. These global tests, rather than the localized Calgary data, support the atlas correctness claim.

Failure behavior is also part of acceptance. Tests deliberately provide insufficient output capacity, require an explicit exception rather than truncation, and then reuse the same architecture instance successfully. Invalid levels, mixed or noncanonical public hierarchical indices, duplicate seeds, empty input, and repeated queries are tested in the same suite. The recorded CUDA memory and race checks complete without reported errors for the earlier GPU release configurations. No equivalent rerun with the four mechanisms enabled is included in the released files.

3.6.7. Separate WGS84 Buffer Experiment

The ellipsoidal extension answers a different timing question from a preindexed fixed hop query. We test the complete 1,482 segment Nose Hill Park boundary at levels 12, 14, and 16 with a 1 km distance. Additional point queries cover levels 20 and 28, the date line, and a near polar location. Each query is run three times in a fresh process using the CPU or GPU halo implementation. Geometric checks use eight OpenMP workers in both variants. The hybrid variant runs on the RTX 4070 with driver 595.79 and CUDA 13.3; the sources are compiled with MSVC 19.51. This boundary buffer study uses a different timing contract from the steady state halo study of Section 4.1.1.

We report reading, source preparation, source covering, implementation initialization, each GPU stage, stopping distance checks, footprint filtering, and output writing separately. Query time includes implementation initialization (GPU initialization in the hybrid variant) but excludes validation and output file writing. CUDA event times are summed over all replay calls. Validation compares complete CPU and GPU rings, and compares the final filtered answer with a spatial hierarchy enumeration that does not use adjacency or the ACM. A separate Python point locator checks 3,000 WGS84 offset samples per query. Additional checks compare 15,840 sampled footprint points with the radius bound, compare pruned and exhaustive source covers in twelve coarse cases, and require capacity exhaustion to fail explicitly. The spatial enumeration shares the distance model, and the Python binding of the PROJ library uses the same geodesic algorithm family; these are complementary checks, not independent proofs of the geodesic solver.

4. Results

We report performance first, followed by component contributions and workload effects, complete workflow costs, and independent topology validation. The preindexed halo study and the WGS84 boundary study have different timing contracts. Their speedups are not interchangeable, and every performance configuration must pass its acceptance checks before its measurements are included.

4.1. Exact Halo Construction

This group measures the reusable ring component on indexed seeds. GPU readiness includes output transfer to pinned host memory; returned container delivery is reported separately. Source covering and geometric distance tests are excluded.

4.1.1. Primary Performance

The first question is whether the GPU improves halo construction when the CPU uses a comparable sparse representation. We therefore compare the GPU implementation with the faster of the serial and twenty thread flat table implementations for each configuration. The primary GPU timer ends when the public output indices reach reusable pinned host memory. The CPU timer includes its fresh returned container, and the GPU performs an additional returned container copy after its timer. These are useful readiness measurements, but not call latencies with equal output ownership. Table 7 reports this comparison by realized output size, and Figure 12(a) shows each configuration.

The GPU is faster in 48 of 81 configurations. The median ratio of the faster flat CPU time to GPU readiness is $2.38$, with range $0.10$ to $9.45$. All 33 CPU wins produce at most 3,000 halo cells. Median ratios rise from $2.46$ for outputs between 3,000 and 40,000 cells to $9.10$ between 0.4 and 1.2 million cells. Dataset medians are similar (Table 8), but this alone does not establish that spatial arrangement has no effect.

The largest query contains 100,000 building seeds at level 28 and eight hops. It returns 10.8 million cells, with selected median GPU readiness of 13.267 ms against 117.132 ms for the twenty thread flat CPU, a ratio of $8.83$. The serial flat CPU takes 715.408 ms; the partitioned and serial container diagnostics take 652.530 and 3529.994 ms. The stronger flat baseline therefore changes the interpretation of acceleration substantially. Table 9 reports the principal configurations without using the slower container reference to headline the GPU result.

The primary summary selects the smallest eligible session median independently for each implementation and configuration (Section 3.6.4); 73 of the 81 CPU/GPU pairs consequently combine different sessions. It describes the best observed eligible performance, not an unbiased estimate of ordinary latency or a simultaneous matched comparison. As a sensitivity check, using each of the two complete matched sessions (A4s1 and A4) instead gives median ratios of $3.13$ and $3.38$, with 49 GPU wins in each. Those all configuration sensitivity summaries retain rows flagged by the spread criterion. They are not an additional clean sample or confidence interval.

Readiness and returned container latency.

The 13.3 ms headline includes the 86.4 MB public output transfer but stops before the query interface copies those indices from the pinned host buffer into its returned vector. The supplied analysis estimates 26.9 ms by adding a mean post timer copy from another session to the selected readiness median. This is a reconstructed estimate, not a measured median call latency. In the mechanism removal session (E2B), readiness and post timer copy have means 13.326 and 14.692 ms, whose sum is 28.018 ms. Table 13 keeps those stages explicit. Neither this sum nor an isolated device stage should be substituted for the primary median without changing the timing interpretation.

CPU thread scaling.

The flat table comparison also permits a direct check of CPU parallelism. On the largest query, the twenty thread time is about $6.1$ times smaller than the serial flat time in the thread sweep. The partitioned implementation remains slower in that session. Conversely, the 1,046 cell coarse query becomes slower when parallel coordination is added. Table 10 and Figure 12(c) show measured scaling; no memory latency or occupancy profiler was run, so these curves are not a diagnosis of the limiting hardware resource.

4.1.2. Component Ablations

The overall ratio combines index conversion, discovery, ring scheduling, and output delivery. We next separate those effects. The dual index study (Section 3.6.2) isolates conversion placement; the schedule study compares broader architectures; and the controlled variant sessions change individual mechanisms or explicitly identified bundles. These studies have different reference implementations and must not be pooled into one component speedup.

In the CPU conversion study with 100,000 seeds and four hops, retaining the compact integer triangular index reduces mean time from 41.753 to 12.198 ms at level 16, a $3.42\times$ ratio, and from 51.205 to 13.164 ms at level 28, a $3.89\times$ ratio. Both variants return 3,124 halo cells in rings of 777, 779, 783, and 785 cells, generate 307,017 logical references, and have equal deterministic hashes. This compares one time boundary conversion with repeated conversion inside our CPU expansion. It is not a measured comparison with BIM.

Alternative execution schedules.

The schedule study of Section 3.6.2 compares host managed sets and sorting schedules against the earlier GPU release of Section 3.4.6. Table 11 retains every ratio rather than summarizing the broad variation by a single mean. These alternatives change several scheduling and set management operations together. Their agreement on cumulative cell sets supports the comparisons, but the ratios do not isolate the necessity of any one component or quantify the improvement over the GPU implementation.

Controlled GPU variants.

The controlled variant sessions measure architectural alternatives and removal of the four mechanisms against their own session references ((b)). Launching a kernel for every ring costs $1.05$ times the reference on the largest query but $2.25$ times on the small coarse query. Replacing warp aggregated append with a per thread counter gives ratios from $1.02$ to $1.16$ across the four listed queries. Omitting the mixing function is particularly poor on the measured coarse case, while block local Morton placement ranges from $0.97$ to $3.45$. These are workload dependent effects, not universal constants.

The 0.7 load target changes readiness by less than five percent on the listed cases, whereas the 0.25 target is up to $1.35$ times slower. Scatter clearing does not improve their readiness appreciably. A 0.7 load target does not by itself imply that all tables can be halved without cost: capacities are discrete, and contention and failure behavior can differ on other workloads. Removing arrival edge skipping has little effect on one listed query and larger effects elsewhere. The preflight and control clearing variant changes two operations together, so their individual contributions remain unresolved.

Output delivery and stage costs.

The largest changes include work outside adjacency. Replacing reusable pinned output with a fresh pageable container costs $2.31$ to $2.76$ times as much for the three large queries in Table 12. This row changes both memory reuse and the delivery boundary. Table 13 and Table 14 report direct stage means, and Figure 12(d) displays them. For the largest query, the traversal stage averages 7.832 ms and output encoding 1.429 ms. Returning a fresh container after pinned readiness adds a measured 14.692 ms in this session.

With all four mechanisms removed, the largest query averages 38.210 ms, including fresh output allocation and delivery. With them, pinned readiness averages 13.326 ms, a ratio of $2.87$. Adding its post timer container copy gives 28.018 ms and reduces the comparison to about $1.36$. The latter is a sum of measured stage means, not a separately timed complete call. Allocation and initialization account for 14.557 ms in the variant, but no page fault counters or GPU idle measurements were collected. The evidence supports reducing output management cost without attributing it to a measured page fault mechanism.

4.1.3. Resolution, Hop Depth, and Output Sensitivity

The number of distinct cells reached, rather than the resolution integer alone, explains much of the change in work. Median selected ratios against the faster flat CPU are $0.33$, $5.19$, and $5.10$ at levels 12, 20, and 28, and $0.47$, $2.91$, and $4.31$ for one, four, and eight hops. These are overlapping descriptive groups, not isolated causal effects of level or depth.

Mapping 100,000 building points produces only 1,914 unique cells at level 12, whereas nearly all remain distinct at levels 20 and 28. Away from exceptional base mesh vertices, an isolated seed contributes $3K(K+1)/2$ cells through K rings after excluding itself. At eight hops this is 108 cells per seed, matching the largest 10.8 million cell output. The realized count is smaller when neighborhoods overlap.

At level 28, 25 of 27 cases attain that isolated count. The 100,000 seed street cases return 2,999,974 cells at four hops and 10,799,646 at eight hops, slightly below 3,000,000 and 10,800,000. At level 20 the median ratio to the isolated count is 0.998 and the minimum is 0.382; level 12 ratios range from 0.008 to 0.362. Thus the fine resolution study mostly measures separate small neighborhoods around sampled points, whereas coarser cases provide substantial overlap. These counts do not measure concurrent contention, and they do not justify treating every level 28 rejection as a return to the parent cell.

4.1.4. Throughput, Memory, and Transfers

Throughput needs a work definition because arrival edge skipping avoids lookup for some geometrically constructed neighbors. We use the logical work W in Equation (20), which counts three neighbor constructions per expanded cell. Dividing by the selected latency yields median rates of 60, 98, and 273 million logical references per second for the serial flat CPU, twenty thread flat CPU, and GPU readiness, respectively. For the 17 cases with at least one million logical references, the medians are 40, 241, and 2,059 million per second. The maximum GPU logical rate is 2,725 million per second.

These are not rates of hash probes or atomic instructions. The largest query has 25.5 million logical references, but arrival skipping suppresses 8.4 million lookup candidates, leaving 17.1 million candidates that enter discovery. Across the primary matrix, counting those lookup candidates instead gives GPU readiness rates with median 189 million and maximum 1,827 million per second. A lookup may itself inspect several slots. Distinguishing these counts avoids attributing skipped work to executed atomic throughput.

The maximum reported GPU workspace is 453,635,624 bytes, or 432.62 MiB. Its reporting formula in Equation (19) contains a $2^{25}$ slot table of 268,435,456 bytes, two output arrays of 86,400,000 bytes each, 10,800,000 arrival bytes, 1,600,000 input bytes, and 168 bytes for reported offsets and control state. A source audit finds that the allocated control structure is 128 rather than 96 bytes, so the explicit device payload is 32 bytes larger than the reported formula. Neither quantity is a measurement of total process or allocator memory. Pinned host staging and the returned vector are additional host allocations.

The largest building query occupies 10.9 million table slots, a load of 0.325. The maximum input and output transfers are 0.763 and 82.398 MiB. The flat CPU implementations report table storage of $8M$ bytes, but do not instrument comparable complete peak memory; we therefore do not claim equal CPU and GPU footprints. Allocation scales with query capacities rather than the more than $4.3\times {10}^{18}$ cells in the level 28 globe. Large distance bands can nevertheless exceed those finite capacities.

4.1.5. Routing and Mapping Costs

The small output cases favor the serial flat CPU, while larger cases provide enough work to amortize GPU readiness costs. A simple local rule chooses the GPU when $r\ge 20$ and the automatic output budget $B_{\mathrm{auto}}=n_{\mathrm{in}}3K(K+1)/2$ of Equation (17), which equals $\left|P\right|3K(K+1)/2$ in this matrix because the input is deduplicated before the experiments, is at least 5,000 cells; otherwise it chooses the serial flat CPU. Available memory and valid capacity bounds remain prerequisites. The level condition summarizes this Calgary sample matrix and is not a mathematical condition for GPU efficiency. The extent occupied by the input cannot cap the output: halos may extend beyond it.

The rule selects the faster implementation in all 81 configurations of the selected session summary. With the same level condition fixed, fitting the threshold on two datasets and testing on the third also selects correctly in all 27 held out cases for each dataset. Several tested thresholds tie, and the stated grid selects 4,000 in these folds. This checks the threshold within a rule family chosen using the same matrix; it is not independent validation of the family, the level condition, a new geography, or different hardware.

Matched session sensitivity is less absolute. The fixed 5,000 rule is correct in 80 of 81 configurations in each complete GPU session. Each discrepancy is a coarse eight hop building query with near equal CPU and GPU times. Thus the perfect selected summary score should not be interpreted as a guarantee. There are no fine level output estimates strictly between 3,000 and 30,000 cells in this matrix, so the precise crossover remains weakly resolved.

Point preprocessing, including assignment to cells, sorting, and deduplication, has median 6.512 ms and maximum 78.205 ms, and exceeds GPU readiness in all 81 corresponding cases. Adding the same recorded preprocessing cost to both implementation times gives a reconstructed median ratio of $1.02$ and maximum $2.16$ against the faster flat CPU. This is not a directly timed paired end to end experiment. It shows little median advantage for one time point workflows while retaining gains in some large cases. Storing indices for repeated queries amortizes that cost. File input, source geometry preparation, and the post timer GPU container copy remain outside this calculation.

Denser point sampling alone does not establish an exact polygon or line cover. The applied matrix evaluates indexed point sets; the separate boundary experiment below evaluates a supplied curve cover and distance predicates.

4.2. Complete Boundary Buffer Queries

This group adds the source geometry and distance tests to the ring component. We first check the spatial interpretation and coverage bounds, then compare complete WGS84 queries. The measurements include geometric processing and implementation initialization; they are not the steady state halo timings of the preceding group.

4.2.1. Spatial Scale and Coverage Checks

We now return from the cost of halo construction to the spatial buffer request that motivated it. A fast ring query is useful here only when its relation to distance and coverage is clear. We first examine sampled adjacent centroid distances, then the base mesh radius bound and the planar illustration. The following subsection reports the separate WGS84 boundary calculation. These forms of evidence answer different questions and must not be treated as interchangeable distance guarantees.

Across buildings, streets, and boundary data, the median level 20 step is approximately 2.495 m. Multiplication by four and eight gives step length totals of 9.98 m and 19.96 m. At level 28, the median step is 0.009747 m and the corresponding totals are 3.90 cm and 7.80 cm. The 5th to 95th percentile interval is 2.24–2.55 m at level 20 and 0.875–0.996 cm at level 28. These are local samples of edge lengths, not exhaustive ring displacement measurements. The exact quantities $L_k$ and $U_k$ are defined in Section 3.2.2 but were not recorded for every performance query.

Accordingly, the four hop building workload is not described as a certified 10 m clearance layer. Dividing 10 m by the level 28 median, which gives about 1,026, is also not a sufficient hop count. Both calculations ignore direction changes and spatial variation. A complete distance query must instead satisfy .

The base mesh audit evaluates the radius constant independently of these samples. Exact rational arithmetic verifies 180 face and vertex radius inequalities and 1,920 supporting plane inequalities for the frozen Pentakis mesh. Its 90 edges have consistent opposite orientations, and all 32 vertex fans are connected. This supports the constant in Equation (27) and the numerical surplus bounds in Table 4. This mesh audit alone does not validate source covering or the numerical distance routine; the separate WGS84 checks are reported in Section 4.2.2.

In the planar point source reference case, the exact point to triangle distance test identifies 364 intersecting cells. The conservative test retains 378 cells and omits none; the first sufficient complete ring is $K_{★}=18$. The farthest retained point is approximately 11.716011 km from the source, below the 11.732051 km cell radius envelope. Across the 24 additional configurations and their 4,205 required cell instances, the check finds no omission or envelope violation. The local boundary illustration in Figure 1 is checked separately: 251 cells intersect its projected 1 km boundary buffer, and the conservative test returns 259 cells with none missing.

4.2.2. WGS84 Boundary Buffer Correctness and Cost

The first stage constructs rings efficiently, but does not alone answer the distance request. We now evaluate both stages together for a WGS84 boundary buffer: complete source covering, adaptive halo calls, whole ring stopping, and cell filtering. All seven query configurations produce identical final sets with CPU and GPU traversal and with the separate spatial hierarchy enumeration. None of the 21,000 offset samples is missing from its returned cell set. All final runs have zero unresolved cells at the stated numerical tolerance. The footprint radius, coarse cover equivalence, and capacity failure checks also pass. Sampling supports the implementation checks but does not replace the geometric bound or certify floating point arithmetic.

Figure 13 shows the actual output of the Nose Hill boundary query at two levels. The finer output follows the requested region more closely. It evaluates the same boundary request as Figure 1, now on WGS84 mapped Pentakis cells rather than the illustrative planar grid. The input is the boundary curve; any central gap is farther than 1 km from that curve and is not a missing filled polygon interior.

Table 15 shows that fast adjacency does not imply equally fast complete distance processing. For the level 16 park boundary, GPU expansion across seven replay calls totals 0.810 ms, while median footprint filtering is 524.63 ms, source covering is 80.83 ms, and GPU initialization is 95.34 ms. Median output writing adds 40.70 ms. The CPU variant is faster for every complete query in this small matrix. Its geometry stage uses the same eight workers, but it avoids GPU setup and transfers. The speedups of the preindexed halo study therefore must not be applied to these boundary buffer times.

The level 16 park query replays depths $1,2,4,8,16,32,64$, totaling 127 ring iterations across seven calls. The median accumulated workspace clearing time is 0.545 ms across those calls, not per call. Accumulated halo call wall time is 3.670 ms; device expansion and hashing account for 0.810 ms. The expansion kernel has three grid wide barriers per hop and one initialization barrier per call, but the recorded stage timer does not isolate barrier cost. Differences between halo wall time and device stage totals include host and runtime work and cannot be attributed to one fixed overhead without profiling.

For Nose Hill at level 16, complete query time is 618.59 ms on CPU and 754.01 ms in the hybrid variant, a $1.22\times$ slowdown. The reported 0.810 ms GPU expansion is only about 0.11% of hybrid wall time, while filtering alone has a 524.63 ms median. Further speeding the adjacency kernel alone therefore cannot materially accelerate this configuration. The 3.69 times slowdown occurs at level 12, where initialization matters more. Point queries are smaller still and pay 75 to 81 ms of GPU initialization for little adjacency work. These fresh process results do not measure a persistent initialized buffer service; subtracting initialization after the fact would not constitute such an experiment. For the measured complete queries, CPU execution is the supported performance choice.

4.3. Independent Topology and Ring Validation

Correctness has two separate questions: whether the atlas describes the intended triangular mesh, and whether a parallel expansion returns the required rings. A shared CPU/GPU atlas cannot independently answer the first. The frozen half edge and refined edge references therefore reconstruct connectivity without using the runtime ACM. The global query tests then concentrate seeds on base mesh edges and around vertices where five or six faces meet. Their scope must also be distinguished from the output agreement checks of the performance harness.

Table 16. Correctness evidence by kind of evidence and implementation scope. The first two reference constructions exclude runtime ACM records. The global ring matrix was run on the earlier GPU release; the performance harness checks cumulative output, not every ring in every timed call.

Evidence	Reference	Coverage	Result
Base transition reference	Frozen oriented half edge representation	81,900 cells at L0–5; 245,700 directed moves; 11,340 directed base mesh edge positions	PASS
Refined edge reference	Canonical refined vertex and edge pairs	245,760 L6 cells; 737,280 directed moves; 11,520 base mesh boundary adjacencies	PASS
Global ring matrix	Exact serial expansion by ring	40 cases; all 180 directed edges; all 32 base mesh vertices; L2,6,12,20,28; $K=1,2,4,8$; 56,160 crossings	PASS
Performance harness output agreement	Exact sorted cumulative CPU/GPU outputs	81 applied configurations with five implementations, plus schedule, mechanism, and thread studies; 1,646 summary rows in total	PASS

Table 16. Correctness evidence by kind of evidence and implementation scope. The first two reference constructions exclude runtime ACM records. The global ring matrix was run on the earlier GPU release; the performance harness checks cumulative output, not every ring in every timed call.

Evidence	Reference	Coverage	Result
Base transition reference	Frozen oriented half edge representation	81,900 cells at L0–5; 245,700 directed moves; 11,340 directed base mesh edge positions	PASS
Refined edge reference	Canonical refined vertex and edge pairs	245,760 L6 cells; 737,280 directed moves; 11,520 base mesh boundary adjacencies	PASS
Global ring matrix	Exact serial expansion by ring	40 cases; all 180 directed edges; all 32 base mesh vertices; L2,6,12,20,28; $K=1,2,4,8$; 56,160 crossings	PASS
Performance harness output agreement	Exact sorted cumulative CPU/GPU outputs	81 applied configurations with five implementations, plus schedule, mechanism, and thread studies; 1,646 summary rows in total	PASS

The global ring matrix passes all 40 configurations with equality of complete cell sets in every ring and deterministic hashes. A separate dense base mesh boundary case with 100,000 seeds and eight hops performs 315,051 logical adjacency operations and 1,447 base mesh edge crossings and returns 5,744 cells. The maximum in the Calgary matrix is 52 crossings, for the level 20 building query with 99,999 unique seeds and eight hops. The largest output case, at level 28, records zero crossings. This characterizes the geographic locality of the applied data rather than validating the atlas.

The ten released measurement files contain 1,646 passing summary rows. This count includes repeated configurations, implementation rows, and controls; it is not 1,646 independent GPU trials. The performance harness compares the final sorted cumulative output after the timed sequence with its CPU reference. Its checks do not verify every ring offset after every timed call. Logical reference counts remain comparable, whereas arrival skipping intentionally changes the number of lookup candidates.

The independent geometry tests still support the unchanged atlas and neighbor equations. However, the released evidence does not show the complete global ring, capacity failure, and device sanitizer suite rerun with all four mechanisms enabled. The insufficient capacity test of the earlier GPU release returns no partial result and permits a subsequent successful call, but does not fill this regression gap. In particular, agreement on cumulative output is not a proof of memory ordering for the cache controlled load described in Section 3.4.6. Per ring tests and an atomic memory access review are required before extending that validation claim to the GPU implementation with the four mechanisms enabled.

5. Discussion

The results answer two connected questions: how efficiently exact rings can be constructed, and what additional work turns that expansion into a spatial buffer query. This section explains the gains in halo construction, the geometric costs that remain, the extent of base mesh portability, and the limits of the available evidence. Exact topological membership, conditional geometric coverage, and numerical implementation checks support different conclusions.

5.1. Why Storage Does Not Grow with the Grid

The first discussion question is which structural choices permit queries at level 28 without global grid storage. The result is not limited to the difference between a GPU and a CPU hash table. Three properties work together. The public hierarchical index supports hierarchy through bit operations; the compact integer triangular index converts local adjacency into arithmetic; and the Atlas of Connectivity Maps isolates base mesh boundary transformations in $3F$ records. The parallel frontier allocates state according to input, hash, and output capacities rather than the complete refined grid. Large requested outputs still require correspondingly large storage.

The geographic inputs have different arrangements, but refinement changes their discrete overlap. At level 28, 25 of 27 configurations reach the isolated seed output count exactly, and the two remaining street cases differ only slightly. These workloads therefore test mostly independent small neighborhoods around geographic samples, rather than connected fine resolution city features. Level 12 and the controlled compact workloads provide overlapping cases; level 20 includes both nearly isolated and substantially overlapping cases. Similar speedups across dataset names alone do not establish independence from spatial morphology.

5.2. GPU Behavior on the Tested Platform

The stage measurements help explain why large and small queries behave differently. For the small cases, initialization, launch, synchronization, encoding, and transfer remain visible even when few cells are returned. A serial flat table avoids those device setup costs. For larger queries, traversal becomes the largest measured device stage, while output encoding and delivery also increase with the number of returned cells (Table 13).

The controlled variants support retaining cooperative ring execution, effective key mixing, and aggregated output insertion for the measured workloads. They also show that reusable pinned staging changes latency substantially, while moving part of returned container delivery outside the readiness timer. Neither the stage means nor the thread sweep measures cache misses, page faults, occupancy, or hardware saturation. We therefore do not interpret the observed ratios as a throughput asymptote or attribute CPU scaling to a demonstrated random access bottleneck. The local routing rule summarizes measured regimes, and requires recalibration before use on different hardware or input arrangements.

5.3. What Complete Distance Queries Add

A fast adjacency operation is necessary for large ring queries, but it is not a substitute for a distance predicate. The WGS84 experiments make this distinction visible: source covering and cell filtering dominate the tested boundary computations, while GPU initialization also matters for small requests. The benefit of the two stage formulation is that the exact halo component can be evaluated and reused independently, while the distance component makes the source geometry, Earth model, stopping rule, and numerical allowances explicit. The measured preindexed readiness gains establish the efficiency of the ring component, not a speedup for every complete spatial buffer. Adding point preprocessing alone reduces the median reconstructed ratio to $1.02$ (Section 4.1.5). The seven WGS84 cases were not rerun with the four mechanisms enabled, so their CPU advantage describes the recorded boundary implementation rather than a comparison with the GPU implementation of the primary matrix.

The measured mapping cost motivates several optimization opportunities that were not evaluated in this study. Batched face plane evaluation could map points in parallel, while coarse spatial preclassification could reduce face searches for nearby points. Fusing mapping, index conversion, and seed insertion could avoid an intermediate host array; asynchronous transfer could overlap preparation of successive batches. For complete boundary queries, persistent GPU initialization, tighter geometric bounds, and further parallel distance processing could reduce the separately measured setup and geometry costs. These are proposed changes, not measured speedups.

5.4. Using a Buffer to Restrict Another Dataset

A complete distance cover can also be used to restrict another spatial dataset. For example, selecting agricultural land near a shoreline requires combining the distance cover with a cover of that land. Let $\mathcal{A}$ be the agricultural land geometry and $A_r$ a cell cover containing every cell that touches it. The cell result is $A_r\cap Q_{r,d}$. Every point of $\mathcal{A}\cap B_d\left(\mathcal{S}\right)$ is covered by this intersection. Traversal must not be restricted to agricultural cells: roads or other land between the shoreline and a field must not block the search. The answer consists of cell candidates for agricultural land. A cell touching a narrow field can still contain a large proportion of nonagricultural land; the geometric guarantee limits its size and distance, not that proportion. Returning entire cadastral parcels is a different query; a selected parcel can extend much farther than $d+\eta$ and is not covered by the surplus bound.

5.5. Portability Across Triangular Base Meshes

The next question is which claims depend on Pentakis and which belong to the execution architecture. The architecture depends on four geometric conditions, not on Pentakis geometry: base mesh faces are oriented triangles, refinement is congruent one into four, every refined cell has three edge neighbors, and each base mesh edge has exactly two incident faces. Given any watertight oriented triangular spherical mesh satisfying these conditions, a preprocessing step enumerates the directed side of each base mesh edge for each incident face and emits $3F$ records $T(f,e)=(g,e^{\prime },\rho )$. Within the supported limits of at most 64 base mesh faces and levels 0 through 28, changing the base mesh changes the atlas but not the current public or compact index layouts, atomic insertion, cooperative frontier execution, or output conversion. Larger face counts or higher levels require a revised packing design, although the local adjacency and atlas principles remain applicable. This establishes a structural adaptation of the method, not a measured implementation or performance result on another base mesh.

Table 17 makes the specialization explicit. The base mesh changes the face selection geometry, projection distortion, number of atlas records, and pattern of boundary crossings. Within these packing limits, it does not change the number of edge adjacent cells for a refined triangle, the local coordinate equations, atomic discovery, frontier scheduling, or output conversion. Pentakis remains our primary experiment because its 60 smaller charts contain vertices where either five or six base mesh faces meet, allowing both arrangements to be tested. Performance portability is not established experimentally in this paper. A Pentakis and icosahedron experiment should match mean cell area, seed geometry, spatial halo reach, and materialized output. Since cell counts are $604^{r_P}$ and $204^{r_I}$, comparable mean area requires $r_I\approx r_P+{\log }_{4}3$, with interpolation between adjacent integer levels or an explicitly matched spatial radius predicate. Equal numeric levels would combine base mesh face count with spatial resolution and would not isolate execution cost.

5.6. Rhombic Macrocharts

The macrochart hypothesis was defined before the ablation matrix: grouping adjacent triangles may reduce between domain transitions without changing the DGGS cells. In the schedule study, the sorting variant with macrocharts has a median cost of $21.24\times$ the earlier GPU release, compared with $20.25\times$ for sorting on natural triangles (Table 11). Thus the tested grouping does not offset sorting and scheduling after each ring. This does not rule out a different macrochart implementation with device resident hashing; that combination was not measured.

5.7. Limitations and Validity

The final discussion question is which conclusions the current experiments do not establish. The following limitations separate exact preindexed halo claims from geometry cover, comparison across DGGSs, and hardware generalization questions.

Input semantics. The building, street, and boundary datasets are converted into deterministic point distributions. Buildings do not contribute every intersected footprint cell, and line sampling is not an exact spherical polygon cover. The halo is exact relative to the indexed cell set, not necessarily to every seed geometry.

Distance completeness. The fixed hop interface alone does not certify a kilometer based query. The extension in Section 3.5.2 integrates complete boundary covering, adaptive halo replay, and curved footprint tests, with separate results in Section 4.2.2. It retains unresolved cells and uses explicit numerical allowances, but lacks a formal outward rounded interval implementation. Consequently, the mathematical guarantee remains conditional on valid distance bounds. The tests do not establish an unconditional certificate, a filled polygon buffer, or the accuracy of the source curve. A finite sample of adjacent distances cannot replace those requirements.

Large distance requests at fine resolution. Small cells reduce the surplus width but increase both the number of rings and the number of returned cells needed for a fixed distance. For scale, one square kilometer alone contains approximately $8.5\times {10}^9$ average level 28 Pentakis cells. The present fixed level array interface is therefore not a demonstrated way to materialize the entire 10 km band around a large source curve at that level. A practical deployment can choose a coarser level from the acceptable surplus, restrict the requested region, or develop a conservative hierarchy with fine cells only near boundaries. These choices must preserve coverage and cannot be replaced by displaying a small sample of the source.

Projection and Earth model. The indexed traversal uses a spherical DGGS with gnomonic charts. The WGS84 extension interprets the longitude and latitude of those spherical points as geodetic coordinates and evaluates distance on the ellipsoid. Neither evaluation studies a new projection or establishes equal area cells. General filled polygon covering and clipping remain separate from the tested boundary curve workflow.

Benchmark scope. The evaluation includes four CPU implementations, a CPU thread sweep, controlled GPU schedules, and thirteen GPU variant rows with explicit session references. The flat CPU table is substantially stronger than the container and partitioned diagnostics. This does not exhaust possible memory layout, vectorization, allocation, or scheduler improvements. H3, S2, and rHEALPix are not timed, and no optimized BIM implementation is benchmarked. Comparisons with those systems must match source semantics, cell scale, spatial reach, output delivery, and timing scope rather than numeric level alone. The absence of such a matched comparison remains an empirical gap.

Timing scope. Primary GPU readiness includes transfer into retained pinned host storage but excludes the subsequent returned vector copy; CPU timings include their fresh result container. Reported readiness ratios therefore do not establish call latency speedups with equal output ownership. Direct stage means and reconstructed completion estimates are labeled separately. Selecting minimum eligible session medians independently for each implementation does not remove timing selection bias; 73 of 81 primary pairs combine sessions. Matched session sensitivity reduces reliance on that selection, but neither repeated calls nor two GPU sessions establish performance across independent machines. Point preprocessing, file input, geometry covering, and rendering are outside the primary interval.

Correctness scope. Independent mesh references validate the atlas without relying on CPU/GPU equality. The global ring matrix additionally tests exact rings across every directed base mesh edge. The performance measurements with the four mechanisms check cumulative output after a timed sequence, not each ring after every call. The global, failure, and sanitizer tests need a documented rerun with all four mechanisms enabled. The cache controlled load, which is not a C++ atomic load, mixed with atomic slot updates also requires a CUDA memory model review; monotone table occupancy and passing output tests are not substitutes for that review.

Capacity exhaustion. The cooperative query does not resize its hash table during a launch. If the before launch estimate is too small, a device status flag terminates the call without exposing a partial result, after which the same instance can retry with a larger bound or select the CPU. The earlier GPU release has a recorded failure and retry test; the current GPU implementation with all four mechanisms enabled needs the regression rerun noted above. Explicit retry can add latency for unexpectedly dense, deep queries. Segmenting a query into bounded epochs or implementing a restartable growable hash could reduce retry cost. Hash growth alone would not remove finite device memory or output array limits. The recorded buffer experiment uses an explicit output capacity of 200,000 cells; no unmeasured maximum supported park resolution is inferred from it.

Hardware generalization. All performance measurements and routing thresholds come from one CPU, GPU, and operating system profile. We make no quantitative speedup, throughput, or routing accuracy claim for another platform. The calibration procedure can be repeated, but the fitted numbers and reported speedups describe the tested hardware. Results on data center GPUs, integrated GPUs, and different CPU memory systems are required to evaluate routing behavior across hardware classes.

Base mesh generalization. The equations and generated $3F$ atlas apply to other oriented triangular spherical meshes, but the current performance matrix uses Pentakis only. A matched area Pentakis and icosahedron benchmark is required to verify that different face counts, base mesh boundary frequencies, and projection locality patterns do not materially change throughput.

Workload regime. Most level 28 geographic samples have effectively disjoint neighborhoods at the tested depths. The coarse matrix, some level 20 cases, and the controlled compact workloads supply overlap evidence. Connected polygon covers, larger distances, and direct contention measurements are needed before generalizing the high resolution speedups to those workloads.

CPU comparison scope. The flat CPU and GPU implementations share compact indices and sparse discovery, but use different memory ownership and synchronization details. The CPU thread sweep reaches about six times the serial speed on two large cases; it does not demonstrate full processor utilization or a memory latency bottleneck. Thread placement and hardware counters were not recorded. Results should not be extrapolated to processors with more cores or different memory systems.

Stage evidence. The stage measurements of the mechanism removal session directly measure representative initialization, traversal, encoding, transfer, and container costs. It is not a decomposition of every primary median. References for mechanism and architectural variants were collected in separate processes, and the fourth mechanism bundles preflight removal with fused control clearing. Page faults and GPU idle fractions were not measured. Claims about those causes would require additional instrumentation.

Data selection. The prepared subsets are deterministic rather than randomized. Requested point subsets are selected by a deterministic stride through the prepared coordinates, and the municipal boundary contains fewer than 100,000 samples. We therefore report actual unique seed counts and avoid interpreting requested point categories as identical work across datasets.

6. Artifacts and Reproducibility Status

The preindexed and WGS84 boundary studies use separate timing contracts, versions, and manifests. The primary C++17 and CUDA matrix uses ten warmups and thirty timed calls per included implementation; its summary draws on four CPU sessions and two GPU implementation sessions. The principal GPU ablation records also use thirty calls, whereas auxiliary CPU controls and the SCALE thread sweep use shorter sequences as specified in the accompanying measurement scripts. The boundary study uses three fresh process repetitions per configuration. Its source covering and distance stages cannot be replaced by the preindexed measurements.

The accompanying material contains the manuscript, editable figures and descriptions, supplied analysis tables, the audited raw measurement summaries, and scripts that check their correspondence. The ten raw files include the E3 CPU baseline and thread studies, ALLOPT and the two A4 matrices, E2B and E2C variants, and earlier E1, E2, and SCALE studies. The numerical audit records the session counter correction, timing ownership distinctions, and the 32 byte workspace reporting discrepancy. Original measurements are retained unchanged. The geometry evidence includes the frozen mesh reference, radius audit, and coverage checks.

The local working source and measurement scripts are available to the authors, but a local manuscript package is not a publicly archived research release. Independent reproduction requires the complete native source, generated atlas and generator, input archives, build and run settings, correctness suite, and raw records. A permanent repository and archival identifier have not yet been assigned. The global regression and memory access checks with the four mechanisms enabled (Section 4.3) also remain prerequisites to a complete release claim. The data availability statement reports this status without inventing a deposit or a formal distance certificate.

7. Conclusions

The principal contribution is an execution architecture for exact halo queries on stored triangular DGGS indices. A public hierarchical index preserves the external representation, while compact integer triangular indices support arithmetic adjacency inside the query. A base mesh atlas stores $3F$ transitions, and a cooperative GPU resident hash frontier retains the cells discovered during expansion. Storage depends on query capacities rather than the complete refined grid. Pentakis is the measured base mesh; performance on other meshes remains untested.

The evaluation uses serial and twenty thread flat table CPU baselines and separates GPU readiness from returned container delivery. On the tested i7 12700K and RTX 4070, the selected session summary gives 48 GPU wins in 81 configurations and a median ratio of $2.38$ against the faster flat CPU. The largest output contains 10.8 million cells and reaches pinned host memory in 13.3 ms, compared with 117 ms for the twenty thread flat table. Its additional returned container copy is not included in that headline. Small queries favor CPU execution; a locally fitted rule describes the current matrix but is not a portable crossover law. Component studies show that discovery, ring control, and output management all affect latency.

Independent mesh references validate Pentakis adjacency. Global tests on the earlier GPU release exercise 180 directed base mesh edges and all 32 base mesh vertices, with 56,160 crossings and exact ring agreement. The performance runs with the four mechanisms agree on cumulative cell sets, but per ring regression and device memory access validation remain distinct requirements. The stronger CPU baseline and explicit delivery scope prevent the GPU ratios from depending on a weak reference or an unstated timing boundary.

A spatial buffer adds geometric obligations that topological speed alone cannot resolve. Complete source covering, sufficient expansion, and conservative cell filtering connect exact rings to distance requests. The separate WGS84 study checks seven configurations, with no missing cells among 21,000 offset samples, but CPU execution is faster in all seven recorded complete workflows. Those cases were not rerun with the four mechanisms enabled. Point preprocessing likewise reduces the reconstructed median advantage to $1.02$. Formal numerical certification, filled polygon semantics, matched spatial comparisons, and further distance processing optimization remain open. The results therefore establish efficient preindexed ring construction and a separate coverage framework, not universally faster complete spatial buffers.