arXiv CSopen access

3DPipe: A Pipelined GPU Framework for Scalable Generalized Spatial Join over Polyhedral Objects

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

databases, sql, data management, storage

3DPipe: A Pipelined GPU Framework for Scalable Generalized Spatial Join over Polyhedral Objects Lyuheng Yuan, Da Yan, Akhlaque Ahmad Indiana University Bloomington {lyyuan,yanda,akahmad}@iu.edu

Fusheng Wang

Stony Brook University [email protected]

arXiv:2604.19982v1 [cs.DB] 21 Apr 2026

Abstract Spatial join is a fundamental operation in spatial databases. With the rapid growth of 3D data in applications such as LiDAR-based object detection and 3D digital pathology, there is an increasing need to support spatial join over 3D datasets. However, existing techniques are largely designed for 2D data, leaving 3D spatial join underexplored and computationally expensive. We present 3DPipe, a pipelined GPU framework for scalable spatial join over polyhedral objects. 3DPipe exploits GPU parallelism across both filtering and refinement stages, incorporates a multi-level pruning strategy for efficient candidate reduction, and employs chunked streaming to handle datasets exceeding GPU memory. Its pipelined execution overlaps CPU data preparation, host-device data transfer, and GPU computation to improve throughput. Experiments show that 3DPipe achieves up to 9.0× speedup over the state-of-theart GPU solution, TDBase, while maintaining excellent scalability. 3DPipe is open-sourced at https://github.com/lyuheng/3dpipe. ACM Reference Format: Lyuheng Yuan, Da Yan, Akhlaque Ahmad and Fusheng Wang. 2018. 3DPipe: A Pipelined GPU Framework for Scalable Generalized Spatial Join over Polyhedral Objects. In Proceedings of Make sure to enter the correct conference title from your rights confirmation email (Conference acronym ’XX). ACM, New York, NY, USA, 15 pages. https://doi.org/XXXXXXX.XXXXXXX

Introduction

Spatial join is a fundamental operation that relates objects across two spatial datasets 𝑅 and 𝑆 based on their geometric relationships. Formally, for each object 𝑟 ∈ 𝑅, spatial join retrieves objects 𝑠 ∈ 𝑆 that satisfy a spatial predicate, including within-𝜏 distance, intersection (i.e., special case with 𝜏 = 0), or 𝑘-nearest neighbors (𝑘-NN). These query types constitute a generalized class of spatial join operations widely used in real spatial applications [19, 38, 44, 47, 51]. Nowadays, many applications require large-scale 3D spatial data processing, including digital pathology [17], human atlases [15], GIS [2], mineral exploration [36], high-definition mapping [49], and urban planning [1]. Meanwhile, advances in AI and sensing technologies are rapidly accelerating the generation of such data. For instance, LiDAR-based autonomous driving systems produce massive Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for components of this work owned by others than the author(s) must be honored. Abstracting with credit is permitted. To copy otherwise, or republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee. Request permissions from [email protected]. Conference acronym ’XX, Woodstock, NY © 2018 Copyright held by the owner/author(s). Publication rights licensed to ACM. ACM ISBN 978-1-4503-XXXX-X/2018/06 https://doi.org/XXXXXXX.XXXXXXX

Figure 1: Multi-Resolution 3D Object Compression point clouds [13, 18], while modern 3D reconstruction [31, 42], generative [45], and segmentation [14, 20, 53] methods enable scalable transformation from raw inputs to structured 3D representations. We adopt polyhedral representations for 3D objects, which are widely adopted in spatial database systems such as PostGIS [6] and PolarDB [5] due to their ability to model complex geometries. Note that point clouds are typically converted to polyhedral representations via surface reconstruction methods [11, 23, 28, 34]. Most works on spatial join focus on 2D spatial objects (minimum bounding rectangles or polygons) [22, 33, 35, 38, 48]. Compared to 2D geometries, 3D objects are substantially more complex, typically consisting of a large number of (triangular) facets, which leads to significant computational overhead. For example, a naïve spatial join approach that directly evaluates object-pair relationships at the facet level incurs prohibitive cost. Such massively parallel facetlevel computations naturally align with the SIMT execution model of GPUs, which is exploited by some recent studies [36, 39, 40]. We are only aware of three works that accelerate spatial join over polyhedral objects: iSPEED [41] which is built with MapReduce, and 3DPro [39] and TDBase [40] which support GPU execution. These works all follow a Filter-and-Refine paradigm, where coarsegrained filtering is first performed using lightweight object-level approximations to remove invalid object pairs, followed by a refinement stage that progressively evaluates candidate pairs using multi-resolution representations derived from geometric simplification (i.e., different levels of detail, or LODs, see Figure 1), enabling early termination whenever the spatial relationship can be determined at lower resolutions, and only falling back to the original high-resolution geometries when necessary. To mitigate the high cost of geometric computations during refinement, 3DPro [39] and TDBase [40] leverage GPU to exploit massive parallelism. While 3DPro [39] and TDBase [40] only use CPU for filtering, we observe that the filtering stage can take up to 82% of the total running time even with OpenMP enabled for multithreaded execution. This is because, besides using object-level minimum bounding box (MBB) for filtering, they also utilize a more precise approximation of object geometry called skeleton-based partitioning, which decomposes each object by grouping nearby facets into intermediate geometric units (i.e., facet clusters) to enable more effective

Conference acronym ’XX, June 03–05, 2018, Woodstock, NY

Figure 2: Illustration of Voxels and Voxel Pairs pruning than object-level MBBs. For convenience, we refer to such facet clusters as voxels, where each voxel is defined as the MBB enclosing the corresponding cluster of facets. We illustrate why skeleton-based partitioning is essential for reducing facet-level computations during refinement using Figure 2, which shows a pair of objects (𝑟, 𝑠), each with two voxels highlighted. To compute the distance between 𝑟 and 𝑠 (i.e., the closest pair of points on their surfaces), denoted by 𝑑 (𝑟, 𝑠), consider a distance 𝑑 1 between a vertex in 𝑟 ’s green voxel and a vertex in 𝑠’s green voxel. Clearly, 𝑑 1 serves as an upper bound of 𝑑 (𝑟, 𝑠). For the orange voxel pair, let 𝑑 2 denote the minimum distance (MINDIST) between their MBBs (efficiently computed by Definition 2 of [37]). Since 𝑑 2 > 𝑑 1 , this voxel pair can be safely pruned, as it cannot contain the closest pair of points that determines 𝑑 (𝑟, 𝑠). By applying this principle across all voxel pairs, a large fraction of voxel pairs can be pruned early, significantly reducing the number of facetlevel comparisons required during refinement. Since voxel-pair filtering introduces much computing overhead, it often becomes a performance bottleneck in TDBase, requiring GPU acceleration. Even in the refinement stage, 3DPro [39] and TDBase [40] do not fully exploit advanced GPU features. For example, they perform distance-bound aggregation in global memory rather than the faster shared memory, and do not leverage CUDA streams to overlap data transfer (of facets) and CPU-GPU computation (of facet-pair distances). In addition, they use a poor implementation that incurs excessive kernel launches, causing substantial scheduling overhead. In this paper, we propose a pipelined GPU framework, called 3DPipe, for scalable generalized spatial join over polyhedral objects. The main contributions are summarized as follows: • 3DPipe exploits GPU parallelism across both voxel-pair filtering and facet-level refinement stages, with (i) chunked streaming to handle datasets exceeding GPU memory, and (ii) a pipelined execution strategy to overlap data transfer, CPU and GPU computations to improve GPU utilization. • 3DPipe adopts a block-centric design, where each GPU (thread) block serves as a basic processing unit to process object pairs (resp. voxel pairs) during voxel-pair filtering (resp. facet-level refinement). Distance-bound aggregations (used for candidate pruning and early termination) are conducted in the fast shared memory using block-wise Hillis-Steele scan [21]. • 3DPipe utilizes a workload flattening strategy to evenly distribute the all voxel-pair (resp. facet-pair) computations for each object pair (resp. voxel pair) among the threads of a block, effectively minimizing the number of kernel launches. • Extensive experiments demonstrate the efficiency and scalability of 3DPipe, achieving up to 9.0× end-to-end speedup over TDBase, the current state-of-the-art GPU solution.

Lyuheng Yuan, Da Yan, Akhlaque Ahmad and Fusheng Wang

The rest of this paper is organized as follows. Section 2 presents the preliminaries, including our offline preprocessing of 3D objects and the used GPU algorithms. Then, Section 3 introduces the detailed design of our 3DPipe framework. Finally, Section 4 reports our experiments, Section 5 reviews the related work, and Section 6 concludes this paper and discusses the future work.

Preliminaries

We explain the preliminaries that are used by our algorithms to be presented in Section 3, including an offline preprocessing stage for object indexing, a brief review of GPU basics and advanced features, and some GPU algorithms used by 3DPipe as basic modules.

2.1

Offline Data Preprocessing

We assume each 3D object is represented as a polyhedron composed of multiple polygonal faces. We also assume each face is a triangle with 3 vertices (called as a facet), since we can always divide a polygon into multiple triangles and regard them as different facets. Given two polyhedra 𝑃1 and 𝑃 2 , we define their distance as 𝑑 (𝑃 1, 𝑃2 ) = min𝑝 1 ∈𝑃1 , 𝑝 2 ∈𝑃2 ∥𝑝 1 − 𝑝 2 ∥ 2 , i.e., the minimum Euclidean distance between any pair of points on their surfaces. Object Voxelization. We preprocess each object by clustering its facets into a set of voxels. Given a target number of voxels 𝑘, we represent each facet by its centroid and apply 𝑘-means clustering over these centroids to group facets into 𝑘 clusters. Each cluster corresponds to a voxel that aggregates spatially proximate facets. For efficient preprocessing, we perform only two iterations of 𝑘-means, i.e., the cluster centroids are updated twice. The initial 𝑘 centroids are uniformly sampled from the vertices of the polyhedron. Object and Voxel Centers. We leverage a simple yet effective upper bound for object- and voxel-level distance-based pruning (recall Figure 2). Specifically, selecting one vertex from each side (either two objects or two voxels) and computing their Euclidean distance yields an upper bound of the minimum distance between the corresponding pair. This bound enables early result confirmation and candidate pruning for object pairs and voxel pairs. To facilitate such computation, we associate each object and each voxel with a representative point, referred to as its anchor (point). For an object, if its MBB center lies inside the polyhedron, the center is directly used as the anchor; otherwise, the anchor is selected as the closest polyhedron vertex to the MBB center. For a voxel, we similarly compute the center of its bounding box, and then select its closest vertex among all facets assigned to the voxel. Given a pair of objects (or voxels), the distance between their anchors provides an efficient upper bound of their minimum distance, which can be computed with negligible overhead (since anchors are precomputed offline) and enables effective pruning. Note that TDBase computes the distance upper bound directly as the distance between bounding box centers, which may lead to incorrect pruning decisions. Figure 3 shows two objects (e.g., blood vessels) with coincident MBB centers while they remain separated, resulting in a non-zero minimum distance but a zero center-tocenter distance. Our anchor selection avoids this issue by ensuring that an anchor always lies on (or within) the object geometry. Level of Detail (LoD). Following 3DPro and TDBase, we adopt a multi-resolution representation by constructing progressively

3DPipe: GPU-Accelerated 3D Spatial Join

Conference acronym ’XX, June 03–05, 2018, Woodstock, NY

blockDim

Object A Object B MBB center 1.0 0.5

Shared Memory

0.0 Z

0.5 1.0

0.5

0.0

0.5

1.0

0.0 Y 0.5 1.0

1.0

Figure 3: Failure Case of Upper-Bounding with MBB Centers simplified polyhedra via iterative mesh simplification [30], where vertices are removed to minimize approximation error and the resulting holes are filled with new facets to preserve connectivity. We materialize a sequence of LoDs along this process, from coarse to fine, with the finest level being the original polyhedron. Simplification strategies differ in geometric guarantees. 3DPro [39] removes only protruding vertices, ensuring each LoD is contained within its higher-resolution counterpart. This monotonic containment yields progressively tighter lower bounds for pruning, but may degrade shape fidelity. In contrast, PPMC [29] minimizes approximation error without enforcing containment, producing higher-quality approximations but lacking valid distance bounds. We follow TDBase [40] and decouple simplification from distance bounding using facet-level Hausdorff bounds (introduced next), enabling progressive tightening of both lower and upper bounds across LoDs while supporting arbitrary simplification algorithms. As in TDBase, we adopt PPMC due to its strong approximation quality. Consistent Voxelization across LoDs. We require voxel partition to remain consistent across all LoDs, since voxel-pair filtering relies on distance bounds computed over voxel pairs (see Figure 2), and its pruning decisions are only valid if the underlying voxelization is invariant throughout the multi-resolution refinement process. We enforce a fixed voxel partition by performing 𝑘-means clustering on the coarsest LoD to determine voxel boundaries, which are then shared across all LoDs. We track the correspondence between facets across LoDs induced by the facet-splitting mesh simplification process, allowing each facet at any LoD to be consistently mapped to its associated voxel. Finally, we compute anchor points for voxels on the original object resolution (highest LoD). Facet-Level Hausdorff Bounds. To enable valid distance bounding under arbitrary LoD simplification, we quantify the geometric deviation between a simplified polyhedron and its original counterpart using Hausdorff distances. Intuitively, given two objects 𝑃1 and 𝑃2 and their simplified versions 𝑃1′ and 𝑃2′ , the true distance 𝑑 (𝑃1, 𝑃2 ) can be bounded by the distance between 𝑃1′ and 𝑃2′ , adjusted by their respective approximation errors. A straightforward approach is to use polyhedron-level Hausdorff distances. However, such bounds are often too loose, as the

Shared Memory

gridDim

Figure 4: GPU and CUDA Concepts

0.5 1.0

Shared Memory

approximation error may vary significantly across different regions of an object (e.g., concave vs. convex areas), leading to ineffective pruning [40]. To obtain tighter bounds, we follow TDBase [40] and refine the error quantification to the facet level. Specifically, for each facet 𝑓 ′ in a low-LoD polyhedron 𝑃 ′ , we precompute two directional distances: (i) Hausdorff distance ℎ𝑑 (𝑓 ′, 𝑃), measuring how far the facet deviates from the original object, and (ii) proxy Hausdorff distance 𝑝ℎ(𝑃, 𝑓 ′ ), which helps derive a tighter lower bound than using ℎ𝑑 (𝑃, 𝑓 ′ ). See [40] for the their detailed definitions. These bidirectional distances allow us to derive tighter lower and upper bounds by considering all facet pairs between two objects. Specifically, [40] proves the following distance bounds: 𝑑 (𝑃1, 𝑃2 ) ≤ 𝑑 (𝑃1, 𝑃2 ) ≥

min

{𝑑 (𝑓1′, 𝑓2′ ) + ℎ𝑑 (𝑓1′, 𝑃1 ) + ℎ𝑑 (𝑓2′, 𝑃2 )},

(1)

min

{𝑑 (𝑓1′, 𝑓2′ ) − 𝑝ℎ(𝑃1, 𝑓1′ ) − 𝑝ℎ(𝑃2, 𝑓2′ )}.

(2)

𝑓1′ ∈𝑃1′ ,𝑓2′ ∈𝑃2′

As the resolution (LoD) of 𝑃 1′ and 𝑃2′ increases, ℎ𝑑 (·) and 𝑝ℎ(·) decrease, yielding progressively tighter bounds on 𝑑 (𝑃 1, 𝑃2 ). In summary, for each facet 𝑓 ′ in each simplified polyhedron 𝑃 ′ (i.e., each LoD), we precompute and store two values offline: 𝑓 ′ .ℎ𝑑 = ℎ𝑑 (𝑓 ′, 𝑃) and 𝑓 ′ .𝑝ℎ = 𝑝ℎ(𝑃, 𝑓 ′ ). These values are later used to compute tight distance bounds during progressive refinement.

2.2

CUDA Concepts and GPU Algorithms

CUDA Basics. The host (CPU) launches a kernel to execute massively parallel computation on the GPU, using the syntax «<gridDim, blockDim»>. As Figure 4 shows, the kernel invocation creates a grid of gridDim thread blocks, where each block contains blockDim threads. Threads are the basic execution units that execute the same kernel function, identified by their intra-block indices threadIdx (starting from 0), while blocks are identified by their indices blockIdx within the grid (also starting from 0). A global thread ID in the grid can be computed as blockIdx * blockDim + threadIdx. On the hardware side, a GPU consists of multiple streaming multiprocessors (SMs), each containing multiple CUDA cores. Thread blocks are scheduled onto SMs for execution. Threads within a block execute on the same SM and can efficiently cooperate via fast on-chip shared memory and barrier synchronization. Execution is further organized in units of warps, each consisting of 32 threads. Threads in a warp execute in an SIMT fashion, i.e., the same instruction is issued to all threads in the warp. A block is thus decomposed into multiple warps, which are scheduled independently by the SM. CUDA Streams. In a typical CUDA program, memory copy between host and GPU (cudaMemcpy) are blocking, resulting in a strictly sequential execution between host-device data transfer and kernel execution. To enable concurrency, CUDA introduces streams, where each stream represents an ordered sequence of

Conference acronym ’XX, June 03–05, 2018, Woodstock, NY

Lyuheng Yuan, Da Yan, Akhlaque Ahmad and Fusheng Wang

stride = 1

stride = 2

stride = 4 INCLUSIVE

–

__shared__ sh_count[1024]; //blockDim = 1024 int tid = threadIdx; for (int stride=1; stride<BlockDim; stride<<=1) { int val = 0; if (tid >= stride) val = sh_count[tid - stride]; __syncthreads(); sh_count[tid] += val; __syncthreads(); } 0 1 3 4 7 9 13 16 EXCLUSIVE

Figure 6: An Illustration of the Hillis-Steele Scan

1: Initialize thread-local aggregator agg 2: for (t ← threadIdx, t < nA·nB, t += blockDim) do i ← t / nB, j ← t % nB 3: Process (ai, bj) to obtain a value val 4: aggregate val to agg 5: Figure 5: Workload Flattening & Thread-Local Aggregation operations (memory copies or kernels) executed on the GPU. Operations within the same stream are executed in issue order, while operations from different streams may overlap. CUDA provides asynchronous primitives such as cudaMemcpyAsync and stream-based kernel launches, allowing data transfer and computation to proceed without blocking the host (which issues the operations to streams). To ensure true asynchrony for memory transfers, host buffers are typically allocated as pinned memory (cudaMallocHost). Dependencies across streams are managed using events. An event can be recorded at a specific point in a stream via cudaEventRecord, which marks the completion of all preceding operations in that stream. Another stream can invoke cudaStreamWaitEvent on this event, ensuring that its subsequent operations will not start until the recorded point is reached. This mechanism enables finegrained cross-stream synchronization without global barriers. Pairwise Workload Flattening. In both filtering and refinement, a thread block is assigned two local sets and needs to process all cross-set pairs between them. This occurs when enumerating all voxel pairs of an object pair in filtering, and all facet pairs of a voxel pair in refinement. Given two sets 𝐴 = {𝑎 0, . . . , 𝑎𝑛𝐴 −1 } and 𝐵 = {𝑏 0, . . . , 𝑏𝑛𝐵 −1 } assigned to a thread block, we need to process all the 𝑛𝐴𝑛𝐵 pairs of their elements (𝑎𝑖 , 𝑏 𝑗 ). To evenly distribute this work among the blockDim threads, we flatten the 2D pair space into a 1D index range in row-major order. Specifically, the pair (𝑎𝑖 , 𝑏 𝑗 ) is mapped to the linear index 𝑡 = 𝑖 · 𝑛𝐵 + 𝑗. As Lines 2–4 in Figure 5 illustrates, the linear indices are assigned to all threads in a round-robin manner to achieve a balanced workload within the block: thread threadIdx processes 𝑡 = threadIdx, threadIdx + blockDim, threadIdx + 2 · blockDim, . . . , until all 𝑛𝐴𝑛𝐵 pairs are covered. For each assigned index 𝑡, the thread recovers the corresponding pair indices by computing 𝑖 = ⌊𝑡/𝑛𝐵 ⌋ and 𝑗 = 𝑡 mod 𝑛𝐵 , and then reads 𝑎𝑖 and 𝑏 𝑗 for computation. Fast Blockwise Aggregation in Shared Memory. Our algorithm frequently needs to aggregate the thread-local results of pairwise computations by taking the minimum or summation. This can be achieved in two steps. In the first step, as shown in Figure 5, we let each thread of a block aggregate all the results of its pairwise computations onto its thread-local aggregator variable agg.

In the second step, the blockDim threads of the block then write their aggregated values from Step 1 onto an array (e.g., sh_count in Figure 6) in the block’s shared memory for fast aggregation. This is achieved by the Hillis-Steele scan as illustrated in Figure 6 for summation aggregation, which completes in log(blockDim) rounds. At round 𝑖, the stride equals 2𝑖 , and each thread tid (when tid ≥ stride) updates its value by adding the element located stride positions before it, thereby maintaining the invariant that, after round 𝑖, each position tid stores the sum over the range [max(0, tid − 2𝑖 + 1), tid] in the initial array. Block-level synchronization barriers are inserted between rounds to ensure that all threads observe a consistent view of shared memory, preventing read-after-write hazards due to independently scheduled warps. After all rounds, each position contains the sum over all preceding elements in the initial array, yielding the final inclusive prefix sum. In particular, the rightmost element stores the final block-level aggregation result. While Figure 6 illustrates the summation case, the same procedure can be directly applied to other associative operations such as minimum aggregation, which is useful for aggregating facet-pair distances (or their bounds) for an object pair to determine its object-pair distance (or its bounds). Notably, TDBase performs such aggregation in global memory which is much slower than shared memory, and multiple threads update a shared result using atomicMin, which leads to significant contention overhead. When the aggregation operator is summation, the inclusive prefix sums produced by the Hillis-Steele scan can be converted to exclusive prefix sums by subtracting the original thread-local value at each position. This variant is particularly useful for computing per-thread output offsets within a block, so that each thread can independently determine its write position in a shared output buffer and avoid write conflicts when emitting multiple results. Triangle-Triangle Distance Computation. We adopt the Möller triangle-triangle distance algorithm [32] to compute the minimum distance between two triangles. The algorithm decomposes the problem into two types of distance computations: (i) point-totriangle distance from each vertex of one triangle to the other triangle, and (ii) segment-to-segment distance for each pair of edges from the two triangles. Concretely, for two triangles, this amounts to six vertex-to-triangle checks (three vertices from each triangle) and nine edge-to-edge checks (3 × 3 edge pairs). The global minimum triangle-triangle distance is guaranteed to be found among these 15 candidate distances [32] by taking their minimum. In the refinement stage, we assign one thread to each facet pair within a voxel pair. Since every thread executes the same Möller routine with the same fixed sequence of 15 candidate distance checks,

3DPipe: GPU-Accelerated 3D Spatial Join

Conference acronym ’XX, June 03–05, 2018, Woodstock, NY

Offline Preprocessing

Object Voxels

Object LoDs, and for each LoD: Object Facets with HD[.], PH[.]

Objects Query (R, S)

MBB-based Object Filtering

vpLB[.], vpUB[.]

Object- and Voxel-Pair Distance Bounding (Algo. 1)

opLB[.], opUB[.]

Object-Pair Pruning

oPairs, r2opOffsets

Voxel-Pair Pruning (Algo. 2)

Voxel-Pair Filtering

Object-Pair Pruning

Result Finalized? NO

YES

Report Result

opLB[.], opUB[.] Facet-Level Distance Bounding (Algo. 4) vPairs

Facet-Level Voxel-Pair Data Preparation

Increase LoD Facet-Level Refinement

Figure 7: Overview of the Execution Pipeline of 3DPipe the computation is highly regular and thus embarrassingly parallel on the GPU. This thread-per-facet-pair design is also adopted by prior GPU-based methods such as 3DPro and TDBase.

The 3DPipe Framework

Our 3DPipe framework supports GPU-accelerated 3D spatial join between 𝑅 and 𝑆 for three types of queries: (1) intersection: which for each object 𝑟 ∈ 𝑅, retrieves all objects 𝑠 ∈ 𝑆 that intersects with 𝑟 , i.e., 𝑑 (𝑟, 𝑠) = 0; (2) within-𝜏 distance: which for each 𝑟 ∈ 𝑅, retrieves all 𝑠 ∈ 𝑆 with 𝑑 (𝑟, 𝑠) ≤ 𝜏; (3) 𝑘-NN: which for each 𝑟 ∈ 𝑅, retrieves 𝑘 objects in 𝑆 that are closest to 𝑟 . In this section, we first overview the execution pipeline of 3DPipe, and then present the details of its voxel-pair filtering, facet-level refinement stages, and additional techniques for 𝑘-NN queries.

3.1

Overview

Figure 7 presents the overall execution pipeline of 3DPipe, where operations involving GPU processing are marked with a GPU icon. • Offline Processing. Each polyhedral object is first preprocessed as described in Section 2.1. This includes (1) voxelizing the object, (2) computing anchor points for the object MBB and its voxels, and (3) progressively simplifying the polyhedron 𝑃 via iterative mesh simplification [29] to obtain a sequence of representations with increasing levels of detail (LoDs), and for each simplified polyhedron (i.e., LoD) 𝑃 ′ , we compute the Hausdorff distance ℎ𝑑 (𝑓 , 𝑃) and proxy Hausdorff distance 𝑝ℎ(𝑃, 𝑓 ) for each facet 𝑓 ∈ 𝑃 ′ . • MBB-based Object Filtering. Given a join query over object sets 𝑅 and 𝑆, we first perform lightweight MBB-based filtering on the CPU using an R-tree built on 𝑆, denoted by 𝑇𝑆 . For a within-𝜏 distance query, for each object 𝑟 ∈ 𝑅, we traverse 𝑇𝑆 from the root and recursively visit a child node 𝑁 only if MINDIST(𝑟, 𝑁 ) ≤ 𝜏. If a leaf node corresponding to an object 𝑠 is reached, we obtain lightweight distance bounds of 𝑑 (𝑟, 𝑠) using their MBBs, denoted by [𝑙𝑏, 𝑢𝑏]. This yields three possible cases: (1) if 𝑢𝑏 ≤ 𝜏, (𝑟, 𝑠) is directly added to the result set; (2) if 𝑙𝑏 > 𝜏, (𝑟, 𝑠) is safely pruned; (3) otherwise, the pair (𝑟, 𝑠) remains undecided. All undecided object pairs form a candidate set oPairs to be subsequently processed by the voxel-pair filtering module. For a 𝑘-NN query, we compute the 𝑘 nearest neighbors for each 𝑟 ∈ 𝑅 over 𝑇𝑆 using a variant of the best-first search algorithm [37]. Specifically, nodes are expanded in ascending order of MINDIST maintained in a priority queue. The search terminates when the

smallest MINDIST in the queue exceeds a threshold 𝜃 . Since only MBBs are considered for objects 𝑠, the exact distance 𝑑 (𝑟, 𝑠) is bounded by an interval [𝑙𝑏, 𝑢𝑏]. Therefore, we set 𝜃 to the 𝑘 th smallest upper bound among the current candidate objects. All candidate object pairs retained in the queue at termination form the set oPairs to be processed by the voxel-pair filtering module. • Voxel-Pair Filtering. Given candidates oPairs, we next perform a finer-grained voxel-level filtering. This stage enables additional object-pair filtering and, for the remaining pairs, applies voxel-pair pruning to reduce the cost of subsequent facet-level refinement. The voxel-pair filtering module has three stage. First, the objectand voxel-pair distance bounding stage utilizes GPU to compute pairwise distance bounds for all voxel pairs of each object pair in oPairs. Specifically, for each object pair, we enumerate all crossobject voxel pairs, compute a lower bound and an upper bound for each voxel pair, and store them in vpLB and vpUB, respectively. At the same time, these voxel-pair bounds are (minimum-)aggregated into object-pair bounds, stored in opLB and opUB. Next, the object-pair pruning stage uses the object-pair bounds opLB and opUB (which are tighter than MBB-based bounds) to prune pairs that can be safely discarded and to directly move those that can be determined to the final result set. Finally, the voxel-pair pruning stage utilizes GPU to prune cross-object voxel pairs whose lower bounds exceed the current object-pair upper bound, since such voxel pairs cannot contribute to the final object-pair minimum distance. The remaining voxel pairs form a reduced set vPairs for subsequent facet-level refinement. • Facet-Level Refinement. Given the reduced voxel-pair set vPairs, we proceed to the final refinement stage to compute exact object distances. Refinement proceeds progressively over LoDs from coarse to fine. At each LoD, the CPU prepares the facet data associated with the current vPairs and transfers them to the GPU. The facetlevel distance bounding stage then computes facet-pair distance bounds (using Eqs. (1)–(2)) for all facet pairs within each voxel pair, and aggregates them to update object-pair bounds opLB and opUB. For 𝑘-NN queries, object-pair pruning is further invoked on GPU to use the tightened object-pair bounds opLB and opUB to prune object pairs that can be safely discarded and to directly move those that can be determined to the final result set. After each LoD, we check whether every object pair can be resolved based on the current object-pair bounds; if so, the results

Conference acronym ’XX, June 03–05, 2018, Woodstock, NY

o2vOffsets

… ℓ𝑥 ℓ𝑦

Lyuheng Yuan, Da Yan, Akhlaque Ahmad and Fusheng Wang

r2opOffsets

ℓ𝑧 𝑢𝑥 𝑢𝑦 𝑢𝑧

op2vpCounts

…

Exclusive Prefix Sum

...

vArray 0

o2vOffsets

oPairs

…

...

vArray 0

(0, 1)

(0, 2)

(0, 8)

(1, 2)

(1, 5)

…

op2vpOffsets

…

Voxel-level Distance Arrays

Object-level Distance Arrays

opLB

vpLB

opUB

vpUB

Figure 8: Data Structures for Object Voxels (Left), Object-Level (Middle) and Voxel-Level (Right) Distance Bounds Algorithm 1 Voxel-Pair Distance Bounding Kernel 1:

(𝑟, 𝑠) ← oPairs[blockIdx]

2: 𝑛𝑟 , 𝑛𝑠 ← voxel counts of 𝑟 and 𝑠 from o2vOffsets 3: offset ← op2vpOffsets[blockIdx] 4: 𝑁 ← op2vpCounts[blockIdx]

//𝑁 = 𝑛𝑟 · 𝑛𝑠

5: for 𝑡 ← threadIdx; 𝑡 < 𝑁 ; 𝑡 += blockDim do

𝑖 ← 𝑡/𝑛𝑠 , 𝑗 ← 𝑡 mod 𝑛𝑠 Load the 𝑖-th voxel of 𝑟 and 𝑗-th voxel of 𝑠 from vArray 8: Compute (𝑙𝑏, 𝑢𝑏) for the voxel pair using their coordinates 9: vpLB[offset + 𝑡] ← 𝑙𝑏, vpUB[offset + 𝑡] ← 𝑢𝑏 10: opLB[blockIdx] ← blockwise minimum of vpLB 11: opUB[blockIdx] ← blockwise minimum of vpUB 6:

are reported immediately. Otherwise, we advance to the next finer LoD and repeat the refinement process until the original-resolution geometry is reached, at which point the bounds become exact object-pair distances and all remaining object pairs can be resolved.

3.2

GPU Algorithms for Voxel-Pair Filtering

Recall that given a query (𝑅, 𝑆), MBB-based object filtering searches R-tree 𝑇𝑆 to identify object pairs (𝑟, 𝑠) that can be directly included in the results. The undecided candidate pairs are then transferred to GPU global memory, where their objects and associated voxels are organized using the data structures shown in Figure 8 (left). Specifically, undecided object pairs are stored in oPairs shown in Figure 8 (middle), with r2opOffsets indexing the candidates associated with each object 𝑟 ∈ 𝑅. The voxels of objects in 𝑅 and 𝑆 are stored in contiguous arrays vArray, with o2vOffsets indicating the starting location for reading voxel coordinates of each object. Each voxel is represented by six coordinates, corresponding to the minimum and maximum values along the 𝑥, 𝑦, and 𝑧 axes. Note that for each object 𝑟 ∈ 𝑅, the number of voxels can be obtained by computing 𝑛𝑟 = (o2vOffsets[𝑟 +1] −o2vOffsets[𝑟 ])/6. We can similarly compute the number of voxels, 𝑛𝑠 , for object 𝑠 ∈ 𝑆. Object- and Voxel-Pair Distance Bounding. For each object pair (𝑟, 𝑠), this stage computes lower and upper distance bounds for all voxel pairs between 𝑟 and 𝑠. The results are stored in vpLB and vpUB, shown in Figure 8 (right). Specifically, the host first computes an array op2vpCounts, which is aligned one-to-one with oPairs. For each (𝑟, 𝑠), the corresponding entry records the number of voxel pairs, i.e., 𝑛𝑟 · 𝑛𝑠 . The host then transfers op2vpCounts to the GPU, where a parallel exclusive prefix sum is performed using an optimized GPU scan primitive (i.e., CUB DeviceScan) to produce op2vpOffsets, which defines, for each object pair (𝑟, 𝑠),

the starting location in vpLB and vpUB for writing the voxel-pair bounds. Given op2vpOffsets, we compute voxel-pair distance bounds on the GPU with gridDim = |oPairs|, as shown in Algorithm 1. Each thread block is assigned to an object pair (𝑟, 𝑠) and enumerates all cross-object voxel pairs via workload flattening. For each voxel pair, we compute its lower and upper distance bounds, stored in vpLB and vpUB, respectively. Meanwhile, these bounds are aggregated within each thread block to obtain object-pair bounds opLB and opUB, using the fast blockwise aggregation technique in shared memory described in Section 2.2. Object-Pair Pruning. Recall that Lines 11–12 have computed tightened object-pair distance bounds [𝑙𝑏, 𝑢𝑏]. For within-𝜏 distance queries, we use CPU to classify the object pairs into three cases: (i) if 𝑢𝑏 ≤ 𝜏, the pair can be directly reported as a result; (ii) if 𝑙𝑏 > 𝜏, the pair is safely pruned; and (iii) otherwise, the pair remains undecided and is forwarded to the subsequent voxel-pair pruning stage for finer-grained processing. For 𝑘-NN queries, object-pair filtering follows a different paradigm since no fixed threshold 𝜏 is given; instead, candidates are maintained and progressively pruned based on relative distance bounds. We defer the discussion to Section 3.4. Voxel-Pair Pruning. Despite prior filtering, the number of voxel pairs within the remaining object pairs can still be substantial. This stage leverages the voxel-pair distance bounds vpLB and vpUB computed by Algorithm 1 to prune voxel pairs that cannot influence the final object-pair distances (recall Figure 2). Specifically, given a voxel pair with distance bound [𝑙𝑏 𝑣 , 𝑢𝑏 𝑣 ] (obtained from vpLB and vpUB) and the aggregated bound [𝑙𝑏𝑜 , 𝑢𝑏𝑜 ] of its corresponding object-pair (obtained from opLB and opUB), we discard the voxel pair if 𝑙𝑏 𝑣 > 𝑢𝑏𝑜 . This voxel-pair pruning stage significantly reduces the workload of the subsequent facet-level refinement, and it is implemented using three GPU kernels which are invoked sequentially. Algorithm 2 shows this three-step procedure for voxel-pair pruning, which is illustrated in Figure 9. Object pairs that have already been pruned in the preceding stage are skipped in this step; their corresponding thread blocks simply return without performing any computation. For simplicity, we omit this detail in Algorithm 2. In the first kernel, each block processes one object pair and scans its voxel pairs to count the number of valid pairs satisfying 𝑙𝑏 𝑣 ≤ 𝑢𝑏𝑜 by workload flattening (Lines 3–5). The per-thread counts are aggregated within the block via shared-memory Hillis-Steele scan to obtain a block-level counter, stored in blkCounter (Line 6). Next, the second kernel performs an exclusive prefix sum over blkCounter to compute blkOffsets using CUB DeviceScan, which determines the starting write position of each thread block in the

3DPipe: GPU-Accelerated 3D Spatial Join

Conference acronym ’XX, June 03–05, 2018, Woodstock, NY

Algorithm 3 Double-Stream Chunked Streaming

oPairs

threshold

opUB

1: Initialize two CUDA streams: compute and memcpy 2: 𝑏𝑒𝑔𝑖𝑛 ← 0, 𝑖𝑡𝑒𝑟 ← 0

(0, 2)

(0, 8)

(1, 2)

(1, 5)

…

op2vpOffsets

…

blkCounter

…

5: while 𝑏𝑒𝑔𝑖𝑛 < |oPairs| do

blkOffsets

…

vpLB 0

globOffset

value

(0, 1)

…

3: Allocate vpLB, vpUB buffers of size CHUNK_SIZE 4: Allocate two output buffers vBuf[0], vBuf[1]

7: 8: 9: 10:

vPairs

…

Figure 9: Voxel-Pair Pruning by Object-Pair Upper Bound

11: 12: 13: 14:

Algorithm 2 Voxel-Pair Pruning 1: // Kernel 1

2: Use blockIdx to obtain (𝑟, 𝑠), 𝑛𝑟 , 𝑛𝑠 , offset and 𝑁 , the same

as in Algorithm 1 Lines 1–4;

𝑢𝑏𝑜 ← opUB[blockIdx]

3: threadCounter ← 0

4: for 𝑡 ← threadIdx; 𝑡 < 𝑁 ; 𝑡 += blockDim do 5:

if vpLB[offset + 𝑡] ≤ 𝑢𝑏𝑜 then

threadCounter ++

6: blkCounter[blockIdx] ← blockwise sum of threadCounter 7: // Kernel 2

8: blkOffsets ← exclusive prefix sum of blkCounter 9: // Kernel 3

10: Perform the same operations as in Lines 2–5

11: localOffsets ← exclusive prefix sum of threadCounter 12: localOffset ← localOffsets[𝑡ℎ𝑟𝑒𝑎𝑑𝐼𝑑𝑥] 13: for 𝑡 ← threadIdx; 𝑡 < 𝑁 ; 𝑡 += blockDim do 14: 15: 16: 17: 18: 19:

𝑙𝑏 𝑣 ← vpLB[offset + 𝑡] if 𝑙𝑏 𝑣 ≤ 𝑢𝑏𝑜 then 𝑖 ← 𝑡/𝑛𝑠 , 𝑗 ← 𝑡 mod 𝑛𝑠 // Get Object IDs globOffset ← blkOffsets[blockIdx] + localOffset vPairs[globOffset] ← (𝑖, 𝑗) localOffset ++

global output array vPairs. Finally, in the third kernel, each thread block re-scans its voxel pairs and writes the surviving pairs to vPairs. Within each block, the valid voxel pairs are first counted by per-thread threadCounter values (Line 10), and the resulting localOffsets are computed via a shared-memory Hillis-Steele scan (Line 11). Each thread then uses its localOffset to determine the write position, and the final global offset is computed as globOffset = blkOffsets[blockIdx] + localOffset. This design ensures that all valid voxel pairs are written into a contiguous array without conflicts, enabling coalesced memory access. Chunked Streaming. The major memory cost of voxel-pair processing comes from the voxel-pair distance-bound arrays vpLB and vpUB, which can be prohibitively large. For example, our experiments show that the total number of voxel pairs can reach 15 billion (see Figure 17), requiring about 15 × 2 × 4 = 120 GB to store vpLB and vpUB, exceeding the memory capacity of commodity GPUs. To address this issue, we process object pairs in a chunked manner. Specifically, we maintain fixed-size buffers on GPU for vpLB and

𝑒𝑛𝑑 ← 𝑏𝑒𝑔𝑖𝑛, 𝑣𝑜𝑥𝐶𝑛𝑡 ← 0 𝑐𝑢𝑟𝑟 ← 𝑖𝑡𝑒𝑟 mod 2, 𝑝𝑟𝑒𝑣 ← 1 − 𝑐𝑢𝑟𝑟 while 𝑒𝑛𝑑 < |oPairs| and 𝑣𝑜𝑥𝐶𝑛𝑡 + op2vpCounts[𝑒𝑛𝑑] ≤ CHUNK_SIZE do 𝑣𝑜𝑥𝐶𝑛𝑡 ← 𝑣𝑜𝑥𝐶𝑛𝑡 + op2vpCounts[𝑒𝑛𝑑] 𝑒𝑛𝑑 ← 𝑒𝑛𝑑 + 1 if 𝑖𝑡𝑒𝑟 > 0 then (Async) D2H copy of vBuf[𝑝𝑟𝑒𝑣] via memcpy stream (Async) launch Algorithms 1 and 2 on compute stream, store results in vBuf[𝑐𝑢𝑟𝑟 ] 𝑏𝑒𝑔𝑖𝑛 ← 𝑒𝑛𝑑, 𝑖𝑡𝑒𝑟 ← 𝑖𝑡𝑒𝑟 + 1

vpUB, and iteratively load a subset of object pairs whose total number of voxel pairs fits within the buffers. For each chunk, we invoke Algorithms 1 and 2 to compute and prune voxel pairs, and send the surviving pairs back to the host to be appended to vPairs. A naïve alternative is to rely on unified memory (allocated via cudaMallocManaged) to transparently migrates data (i.e., vpLB, vpUB, vPairs) between host and device memory. However, unified memory operates at the granularity of pages, whereas voxel-pair processing exhibits large, sequential access over GB-scale arrays. As a result, data is migrated in a page-by-page manner, leading to a massive number of page faults that incur severe overhead. Double-Stream Pipelining. While chunked streaming effectively bounds GPU memory usage, a naïve implementation suffers due to synchronous device-to-host (D2H). In particular, after processing each chunk, the valid voxel pairs must be copied back to host memory before the next chunk can be processed, leaving the GPU idle during the transfer. To further improve throughput, we overlap D2H transfer with GPU computation using two CUDA streams: a compute stream for kernel execution and a memcpy stream for asynchronous output transfer. Algorithm 3 provides the details. Specifically, Line 3 allocates a single pair of vpLB and vpUB buffers for use by the compute stream in Line 13, which is sufficient since the compute stream processes one chunk at a time. Note that Lines 6 & 8–10 construct the next chunk of object pairs from oPairs so that the total number of their voxel pairs does not exceed CHUNK_SIZE, hence vpLB and vpUB buffers will not overflow. Line 4 allocates two output buffers vBuf[0] and vBuf[1] to store the valid voxel pairs from consecutive chunking iterations, so that when one buffer is used for computation, the other is being transferred. Specifically, at each iteration, we alternate between the two buffers, where curr stores the output of the current chunk and prev holds the output from the previous chunk. When iter > 0, the results in vBuf[prev] are asynchronously copied back to host memory via the memcpy stream (Lines 11–12), while the compute stream simultaneously processes the current chunk and writes its results to vBuf[curr] (Line 13). Since the host issues both operations without blocking, and CUDA allows concurrent execution across streams,

Conference acronym ’XX, June 03–05, 2018, Woodstock, NY compute stream memcpy stream

C1-A1

C1-A2

C2-A1 C1-D2H

C2-A2

Lyuheng Yuan, Da Yan, Akhlaque Ahmad and Fusheng Wang C3-A1

C3-A2

C2-D2H

Figure 10: Pipelined Computation and Memcpy the D2H transfer of one chunk overlaps with the computation of the next chunk. This design improves GPU utilization. The double-buffering scheme is implemented using two events to coordinate the two CUDA streams. First, after the compute stream finishes producing vBuf[curr] for the current chunk in Line 13, it records a computeDone event. Before the memcpy stream issues the D2H transfer for that buffer in Line 12, it waits on this event to ensure that the output has been fully generated. Second, after the D2H copy of vBuf[prev] completes in Line 12, the memcpy stream records a copyDone event. When the same buffer is reused by a later iteration, the compute stream waits on the corresponding copyDone event before writing new results into it in Line 13, thereby avoiding overwrite of data that is still in transit. Figure 10 illustrates the pipelined execution of computation and D2H memory transfer. Here, “C1” denotes chunk 1, and “A1” (resp. “A2”) corresponds to Algorithm 1 (resp. Algorithm 2). Summary. Although the filtering stage accounts for only a small fraction (e.g., 15%) of the total runtime after GPU acceleration, it is critical for enabling low-latency coarse search. In practice, it produces the majority of final results (up to 90%) within a few seconds, while refining the remaining undecided object pairs is significantly more expensive. This makes our GPU-based filtering particularly valuable in scenarios requiring partial results ASAP.

3.3

Facet-Level Refinement with GPU

The filtering stage in Section 3.2 produces valid voxel pairs which are further evaluated in the refinement stage. Recall from Figure 2 that each voxel contains multiple triangle facets in close proximity. The refinement stage computes pairwise facet-level distances for all valid voxel pairs to obtain precise object distances. Recent GPU solutions 3DPro [39] and TDBase [40] adopt the Möller algorithm [32] for triangle-triangle distance computation as described in Section 2.2, to exploit its regular structure and massive parallelism. Specifically, they assign each triangle pair to a thread, and every thread executes the same fixed sequence of 15 candidate distance checks. This design leads to an SIMT-friendly computation. Despite this parallel formulation, we identify two performance bottlenecks that significantly limit their efficiency. The first bottleneck is excessive kernel launches. Given two voxels 𝑀 and 𝑁 , existing approaches such as TDBase launch one kernel per triangle in 𝑀, where each launch computes the distances between a single triangle in 𝑀 and all triangles in 𝑁 . As a result, computing all triangle-level distances between 𝑀 and 𝑁 requires |𝑀 | kernel launches. This design is inefficient, as frequent kernel launches incur non-trivial scheduling overhead on the GPU. The other bottleneck is the inefficient use of memory during aggregation. After computing all pairwise facet-level distances within a voxel pair, the results must be aggregated to obtain the minimum distance for the a voxel pair. This operation is memory-intensive and ideally suited for GPU shared memory due to its low latency and high bandwidth. However, existing approaches such as TDBase

Algorithm 4 Facet-Level Refinement Input: vPairs, facet-level Hausdorff distance arrays HD_r, HD_s, and facet-level proxy Hausdorff distance arrays PH_r, PH_s. Output: Voxel-pair distance bounds vpLB, vpUB 1: (𝑣 𝑟 , 𝑣 𝑠 ) ← vPairs[blockIdx] 2: 𝑛𝑟 ← number of triangles in voxel 𝑣 𝑟 3: 𝑛𝑠 ← number of triangles in voxel 𝑣 𝑠 4: localLB ← FLT_MAX, localUB ← FLT_MAX 5: for 𝑡 ← threadIdx; 𝑡 < 𝑛𝑟 · 𝑛𝑠 ; 𝑡 += blockDim do 6: 𝑖 ← 𝑡/𝑛𝑠 , 𝑗 ← 𝑡 mod 𝑛𝑠 7: Read coordinates of the 𝑖-th facet in 𝑣𝑟 and 𝑗-th facet in 𝑣𝑠 8: Compute their triangle-triangle distance 𝑑𝑖𝑠𝑡 9: lb ← 𝑑𝑖𝑠𝑡 − PH_r[𝑣𝑟 ] [𝑖] − PH_s[𝑣𝑠 ] [ 𝑗] // Eq. (2) 10: ub ← 𝑑𝑖𝑠𝑡 + HD_r[𝑣𝑟 ] [𝑖] + HD_s[𝑣𝑠 ] [ 𝑗] // Eq. (1) 11: localLB ← min(localLB, lb) 12: localUB ← min(localUB, ub) 13: Aggregate thread-local bounds in shared memory to obtain block-wise bounds (vPairLB, vPairUB) 14: vpLB[blockIdx] ← vPairLB, vpUB[blockIdx] ← vPairUB perform aggregation directly in the slower global memory by letting threads update current aggregated value by atomicMin, which incurs higher access latency and contention that forces serialization. Facet-Level Refinement. We design an algorithm that launches only a single kernel for efficient facet-pair distance computation and aggregation at the voxel-pair level. Algorithm 4 shows this kernel where each thread block processes one voxel pair and computes its distance bounds based on all its facet pairs. The resulting voxelpair bounds (vpLB and vpUB) are then transferred back to the host, which aggregates them into the object-pair bounds (opLB and opUB) for determining whether each object pair can be pruned or finalized. This refinement is then repeated for the next LoD if needed. In Lines 1–3 of Algorithm 4, each thread block is assigned a voxel pair (𝑣𝑟 , 𝑣𝑠 ) and retrieves the numbers of triangles in the two voxels. The loop at Line 5 enumerates all facet pairs between 𝑣𝑟 and 𝑣𝑠 . Specifically, Line 6 flattens the 2D facet-pair space into a 1D index, so that each thread processes facet pairs in a round-robin manner. For each assigned pair, we read the triangle coordinates (Line 7) and compute their triangle-triangle distance (Line 8). Using the precomputed facet-level Hausdorff and proxy Hausdorff distances, Lines 9 and 10 derive the lower and upper bounds according to Equations (2) and (1). Each thread maintains its local bounds during the iteration (see Lines 11–12). After all facet pairs are processed, thread-local results are aggregated in shared memory to obtain blockwise voxel-pair bounds, which are written to vpLB and vpUB. Facet-Level Voxel-Pair Data Preparation. Algorithm 4 requires accessing several pieces of facet-level data for each valid voxel pair. Specifically, in Line 1, each thread block retrieves a voxel pair (𝑣𝑟 , 𝑣𝑠 ) from vPairs, which encodes the offsets and lengths of the corresponding facet arrays. Lines 2–3 obtain the numbers of triangles in the two voxels, determining the total number of facet pairs to be processed. Inside the main loop from Line 5, each thread reads the triangle coordinates of its assigned facet pair (Line 7) for distance computation. In addition, Lines 9–10 access the precomputed facet-level Hausdorff and proxy Hausdorff values (HD and PH) associated with each facet to derive tight lower and upper bounds.

3DPipe: GPU-Accelerated 3D Spatial Join

Coordinates 0

Hausdorff Distance

(x1, y1, z1, x2, y2, z2, x3, y3, z3) A

vPairs

Conference acronym ’XX, June 03–05, 2018, Woodstock, NY

I 5

…

2: Initialize two CUDA streams

……

3: for 𝑐 𝑖 in chunks [𝑐 0 , 𝑐 1 , 𝑐 2 , · · · , 𝑐 𝑛 ] do

……

Read by Block 1

……

(for current LoD)

Read by Block 0

Format: (vr_offset, vr_len, vs_offset, vs_len)

(x1, y1, z1) (x2, y2, z2)

A (x3, y3, z3)

……

Algorithm 5 CPU-GPU Pipelining in Refinement Stage 1: Launch a dedicated CPU thread for data preparation

(hdJ, phJ)

……

𝑐𝑢𝑟𝑟 ← 𝑖 mod 2, 𝑝𝑟𝑒𝑣 ← 1 − 𝑐𝑢𝑟𝑟 Wait until the facet-level data of 𝑐𝑖 is prepared 6: (Async) perform H2D facet-level data transfer for 𝑐𝑖 , and launch Algorithm 4 for 𝑐𝑖 , both on stream 𝑐𝑢𝑟𝑟 7: if 𝑖 > 0 then Wait for stream 𝑝𝑟𝑒𝑣 to complete the computation on 𝑐𝑖 −1 8: 9: Perform D2H transfer of vpLB and vpUB for 𝑐𝑖 −1 10: Aggregate them to object-pair bounds (opLB, opUB) for 𝑐𝑖 −1 11: Perform D2H bounds transfer and aggregation for 𝑐 𝑛 4: 5:

𝑂1

𝑂2

𝑂3

𝑂4

Filtering

[1.2, 2.5]

[1.8, 3.2]

[2.2, 2.8]

[1.5, 1.9]

Refinement 20% LoD

[1.2, 2.4]

[1.8, 2.5]

[2.2, 2.7]

[1.5, 1.7]

Refinement 40% LoD

[1.2, 2.3]

[1.8, 2.1]

[2.2, 2.5]

Refinement 60% LoD

[1.2, 1.7]

[1.8, 1.9]

Figure 11: Voxel and Voxel-Pair Data Layout H2D

CPU Stream A Stream B

GPU Computation

Data Preparation C2

C5 C3

Figure 13: Progressive Candidate Pruning for a 𝑘-NN Query Figure 12: CPU-GPU Pipelining in Refinement Stage To support efficient GPU execution, these data for the current LoD are first prepared on the CPU and organized into contiguous arrays following the layout illustrated in Figure 11, and then transferred to GPU global memory. Specifically, the CPU first collects all valid voxel pairs produced by the filtering stage, deduplicates their voxels to get a list of unique voxels, and obtains their facets which are then organized into two compact facet arrays, one for facet coordinates and the other for facet-level Hausdorff distances. Each voxel pair (𝑣𝑟 , 𝑣𝑠 ) is then represented by a tuple (𝑣𝑟off, 𝑣𝑟len, 𝑣𝑠off, 𝑣𝑠len ), which records the offsets and lengths of their corresponding facet segments in the two facet arrays. This organization allows each thread block of Algorithm 4 to directly access the required data via simple pointer arithmetic. Chunked Pipelining by CUDA Streams. We observe that the facet-level data preparation phase incurs a non-negligible time cost (nearly 59% of the total refinement time) when data preparation and GPU-accelerated facet-level refinement are executed in sequence. To overlap data preparation with GPU computation, we partition the valid voxel pairs into multiple chunks and process them in a pipelined manner. We implement this pipelining via two independent CUDA streams. Given a sequence of chunks [𝑐 0, 𝑐 1, 𝑐 2, · · · ], the host assigns even-indexed chunks to one stream and odd-indexed chunks to the other stream. The benefit is when one stream uses GPU for computation, another stream may copy prepared data from CPU to GPU, while CPU simultaneously prepares data for a subsequent chunk. As a result, GPU-based computation, CPU-based data preparation, and host-to-device (H2D) data transfer are effectively overlapped to minimize resource waste, as Figure 12 illustrates. Algorithm 5 describes the implementation of the pipelined execution. First, the main thread launches a CPU thread dedicated for preparing data (Line 1), and initializes two CUDA streams (Line 2).

Then, we iterate every chunk 𝑐𝑖 , identify which stream 𝑐𝑖 should use (Line 4). Once the data of 𝑐𝑖 is prepared by the CPU thread (notified via a condition variable, Line 5), the main thread launches H2D facet-level data transfer and then the refinement kernel for 𝑐𝑖 on the current GPU stream (Line 6). Meanwhile, we proceed to receive vpLB and vpUB for 𝑐𝑖 −1 (Line 9) after being computed by the other stream (Line 8, waiting via cudaStreamSynchronize) and aggregate them to object-pair bounds (opLB, opUB) (Line 10). Finally, Line 11 completes this host-side processing for the last chunk. Chunked pipelining also reduces GPU memory consumption: instead of materializing facet-level data for all valid voxel pairs at once, the GPU maintains only the data for two consecutive chunks at any time, one for each stream under a double-buffering scheme.

3.4

GPU Algorithms for Object-Pair Pruning

Recall from Figure 7 that Object-Pair Pruning is used in both stages of Voxel-Pair Filtering and Facet-Level Refinement, which classifies the object pairs into three status CONFIRMED, REMOVED or UNDECIDED based on their distance bounds in opLB and opUB. In Section 3.2, we have seen how to perform Object-Pair Pruning using CPU for within-𝜏 distance queries based on the fixed threshold 𝜏. This section designs a GPU kernel for Object-Pair Pruning in 𝑘-NN queries, which determines the top-𝑘 results progressively as object-pair bounds become tighter across refinement levels, as Figure 13 illustrates. Specifically, consider four candidates 𝑂 1, 𝑂 2, 𝑂 3, 𝑂 4 for a top-2 query. After filtering, the bounds are still too loose to finalize any result: although 𝑂 4 is already known to be closer than 𝑂 3 , its relative order with respect to 𝑂 1 and 𝑂 2 remains unresolved. At 20% LoD, refinement tightens the bounds enough to prove that 𝑂 4 is closer than both 𝑂 2 and 𝑂 3 . Hence, 𝑂 4 must belong to the top-2 result set and is marked as CONFIRMED. At 40% LoD, the updated bounds show that 𝑂 3 is necessarily farther than 𝑂 2 , so

Conference acronym ’XX, June 03–05, 2018, Woodstock, NY

Lyuheng Yuan, Da Yan, Akhlaque Ahmad and Fusheng Wang

Algorithm 6 𝑘-NN Object-Pair Pruning

Table 1: Dataset Parameters Dataset

1: 𝑟 ← blockIdx 2: kLeft ← 𝑘 − numConfirmed[𝑟 ] 3: offset ← r2opOffsets[𝑟 ], 𝑁 ← r2opOffsets[𝑟 + 1] − offset

Digital Pathology

4: for 𝑡 ← threadIdx; 𝑡 < 𝑁 ; 𝑡+ = blockDim do 5: 6: 7: 8: 9:

10: 11: 12: 13: 14: 15:

_ , 𝑚 ← oPairs[offset + 𝑡] if status[𝑚] ≠ UNDECIDED then continue 𝑐𝑙𝑜𝑠𝑒𝑟 ← 0, farther ← 0 for 𝑗 ← 0; 𝑗 < 𝑁 ; 𝑗 + + do _ , 𝑛 ← oPairs[offset + 𝑗] if 𝑛 = 𝑚 or status[𝑛] ≠ UNDECIDED then continue if opLB[𝑚] ≥ opUB[𝑛] then 𝑐𝑙𝑜𝑠𝑒𝑟 ++ if opUB[𝑚] ≤ opLB[𝑛] then farther++ if 𝑁 − farther < kLeft then status[𝑚] ← CONFIRMED else if 𝑐𝑙𝑜𝑠𝑒𝑟 ≥ kLeft then status[𝑚] ← REMOVED else status[𝑚] ← UNDECIDED

𝑂 3 is marked as REMOVED. At 60% LoD, 𝑂 1 is further shown to be closer than 𝑂 2 , which makes 𝑂 1 the remaining top-1 candidate and therefore CONFIRMED. This example shows that 𝑘-NN evaluation can progressively confirm results and discard impossible candidates without waiting for exact distances of all object pairs. While Figure 13 illustrates the process of Object-Pair Pruning for four candidate objects, in large datasets CPU-based top-𝑘 selection becomes a bottleneck (e.g., tens of seconds); we therefore offload it to GPU, reducing the per-round cost to subseconds. Algorithm 6 shows our GPU kernel. We assign one query object 𝑟 ∈ 𝑅 to each thread block, and let the threads in the block evaluate the candidate object pairs associated with 𝑟 . Specifically, Line 1 identifies the query object 𝑟 by blockIdx. Line 2 computes kLeft, the number of result objects that still need to be determined for 𝑟 , by subtracting numConfirmed[𝑟 ] from 𝑘. Here, numConfirmed[𝑟 ] records how many candidates of 𝑟 have already been confirmed as top-𝑘 results in previous LoD pruning rounds (initialized as 0). Line 3 then obtains the contiguous subarray of candidate object pairs of 𝑟 in oPairs using r2opOffsets (recall Figure 8). The for-loop from Line 4 lets the threads examine the candidates of 𝑟 in parallel. Each thread takes one candidate object pair 𝑚 at a time (Line 5), skipping those whose status has already been determined (Line 6). For each remaining candidate 𝑚, we compare its current distance bounds with every other undecided candidate 𝑛 of the same query object 𝑟 (Lines 8–12), where we count how many candidates are guaranteed to be closer (resp. farther) than 𝑚, which is recorded by local variable 𝑐𝑙𝑜𝑠𝑒𝑟 (resp. farther). Lines 13–15 then determine the status of 𝑚 based on these two counters. This pruning procedure is invoked after each update of objectpair bounds, so more candidates can be progressively confirmed or removed as refinement proceeds from coarse to fine LoDs.

Experiments

Besides our 3DPipe framework, only 3DPro [39] and TDBase [40] support 3D spatial join on GPU. However, 3DPro is not opensourced, while TDBase is reported to outperform 3DPro by up to 4× in time [40] and is open-sourced [7]. In this section, we compare 3DPipe with TDBase as the state-of-the-art baseline to evaluate efficiency and scalability. We also conduct ablation studies to assess

ModelNet40

Join Task Nuclei, Vessel Nuclei, Nuclei Train, Test Train, Train

Alias NV NN TI TT

𝑘 [1, 10] [1, 10] [1, 10] [5, 50]

𝜏 [0, 400] [0, 200] [0, 20] [0, 10]

the effectiveness of our optimizations for filtering and refinement. 3DPipe is open-sourced at https://github.com/lyuheng/3dpipe.

4.1

Experiment Setup

Datasets. Following existing works [39–41], we use two real-world 3D digital pathology objects which are generated by reconstructing 2D segmented objects from brain tissues [25, 26], represented in Object File Format (OFF) format [4]. Specifically, one object is a blood vessel containing 30,000 facets and 5 bifurcations; another object is a nucleus cell which has 300 facets. We follow the procedure of [39, 40] to replicate the vessel into multiple objects and shift them to different locations (without overlap of their bounding boxes) to simulate blood vessels in human body. We denote the resulting vessel object dataset as V. We also replicate the nucleus cell and uniformly distribute the cells in the space of V to form a nuclei cell dataset, denoted as N. Unless otherwise specified, we assume |V| = 5000 and |N| = 1 million. We can change spatial join workload by varying the numbers of replicated objects in V and N. We also use ModelNet40 [3, 46], which contains 12,311 prealigned 3D CAD models in OFF format from 40 categories such as beds and chairs. It is split into 9,843 training objects and 2,468 test objects. By replicating each object 100 times and shifting the replicas to different locations, we obtain a training dataset T of 984,300 objects and a test dataset I of 246,800 objects. Queries. We evaluate all three query types: intersection, within-𝜏 distance, and 𝑘-NN, with parameter settings summarized in Table 1. We choose the ranges of 𝑘 and 𝜏 according to the characteristics of each dataset so as to generate representative join workloads without making the running time excessively long; note that 𝜏 = 0 is a special case corresponding to intersection queries. Configurations. All experiments were conducted on a Polaris node at Argonne National Laboratory with 4 AMD EPYC 7543P CPUs (32 cores), 512 GB DRAM, and one NVIDIA A100 GPU with 40 GB memory connected via PCIe Gen4. We use CUDA 12.8 and compile all implementations with -O3. CPU-side processing in both 3DPipe and TDBase uses OpenMP with 32 threads whenever applicable, e.g., to parallelize MBB-based Object Filtering over all objects 𝑟 ∈ 𝑅 against the static R-tree 𝑇𝑆 . For voxelization, the number of voxels is set to 2% of the number of facets per object. For pipelining, the chunk sizes in Algorithms 3 and 5 are set to 5 GB and 500,000 voxel pairs, respectively.

4.2

Comparison with TDBase

Figure 14 compares 3DPipe with TDBase on all three query types, including within-𝜏 distance, intersection (𝜏 = 0), and 𝑘-NN. Overall, 3DPipe consistently outperforms TDBase across all datasets and query settings, with especially large gains on 𝑘-NN queries.

3DPipe: GPU-Accelerated 3D Spatial Join

Conference acronym ’XX, June 03–05, 2018, Woodstock, NY

3DPipe

100

0 50 100 150 200 250 300 350 400 τ

(b) 𝑘-NN NV (varying 𝑘)

(a) 𝑘-NN NV (varying 𝑘)

(b) 𝑘-NN NN (varying 𝑘)

200

0 1 2 3 4 5 6 7 8 9 10 k

Runtime (sec)

0 1 2 3 4 5 6 7 8 9 10

(e) Within-𝜏 TI (varying 𝜏)

Runtime (sec)

(f) 𝑘-NN TI (varying 𝑘) 200

150 100

100

50 0 1 2 3 4 5 6 7 8 9 10 τ

(g) Within-𝜏 TT (varying 𝜏)

(d) 𝑘-NN TT (varying 𝑘)

On TT, shown in Figures 14g and 14h, the speedups further increase to 5.4× for within-𝜏 and 9.0× for 𝑘-NN. These results show that 3DPipe consistently outperforms TDBase across datasets and query types, with larger gains under heavier workloads, confirming the importance of 3DPipe’s fully GPU-oriented optimization.

4.3

200

0 5 10 15 20 25 30 35 40 45 50

Figure 15: Runtime of Filtering Stage

100

0 2 4 6 8 10 12 14 16 18 20

200

0 5 10 15 20 25 30 35 40 45 50

(d) 𝑘-NN NN (varying 𝑘)

300

150

0 1 2 3 4 5 6 7 8 9 10 60

100

0 20 40 60 80 100 120 140 160 180 200 τ

Runtime (sec)

100

0 1 2 3 4 5 6 7 8 9 10

100

0 1 2 3 4 5 6 7 8 9 10

0 1 2 3 4 5 300 6 7 8 9 10

200

0 1 2 3 4 5 6 7 8 9 10

(a) Within-𝜏 NV (varying 𝜏) 300

300

Runtime (sec)

200

Runtime (sec)

200

TDBase 40 Runtime (sec)

150

Runtime (sec)

300 Runtime (sec)

Runtime (sec)

TDBase 300

Runtime (sec)

3DPipe 100

0 5 10 15 20 25 30 35 40 45 50

(h) 𝑘-NN TT (varying 𝑘)

Figure 14: Performance Comparison with TDBase On NV, Figures 14a and 14b show that 3DPipe achieves up to 2.7× speedup for within-𝜏 queries and 5.2×–7.8× speedup for 𝑘-NN queries. The improvement comes from accelerating both filtering and refinement on GPU, rather than leaving key bottlenecks on CPU (TDBase uses only CPU for filtering). The advantage is particularly pronounced for 𝑘-NN, where TDBase spends substantial CPU time on object-pair filtering, while our GPU-based design in Section 3.4 keeps each round of 𝑘-NN processing at the sub-second level. A similar trend is observed on NN in Figures 14c and 14d. 3DPipe reaches up to 3.4× speedup for within-𝜏 queries and 3.9× for 𝑘-NN queries. In addition, our runtime curves are noticeably more stable than those of TDBase. For example, in Figure 14c, when 𝜏 increases to 80, TDBase rises from 22 s to 109 s, whereas 3DPipe increases only from 13 s to 33 s. A similar jump appears in Figure 14d when 𝑘 increases from 2 to 3: TDBase grows from 53 s to 243 s, while 3DPipe increases from 17 s to 67 s. This indicates that 3DPipe scales more robustly as the number of voxel or facet comparisons grows. The same pattern holds on ModelNet40. As shown in Figures 14e and 14f, 3DPipe achieves up to 4.6× and 6.1× speedup, respectively.

Performance Breakdown

Filtering Stage. Figure 15 reports the runtime of the filtering stage for 𝑘-NN queries on NV, NN, TI, and TT datasets. We focus on 𝑘-NN queries, as filtering is a major performance bottleneck in this setting, whereas for within-𝜏 queries the filtering stage completes only within a few seconds and does not dominate the overall runtime. Our filtering stage consists of two phases: (i) MBB-based object filtering on CPU, and (ii) voxel-pair filtering with GPU. Overall, 3DPipe significantly outperforms TDBase, achieving 15×–18× speedup on NV, 4×–6× on NN, and up to 5.3×–6.6× on ModelNet40. The sources of improvement differ across datasets. On NV, MBBbased filtering is lightweight, and the runtime is dominated by voxel-pair filtering. In this case, our GPU-based voxel-pair distance computation and pruning (Algorithms 1 and 2) effectively exploit massive parallelism, while TDBase relies on CPU parallelism (OpenMP), leading to a substantial performance gap. In contrast, on NN, TI, and TT, MBB-based filtering becomes the bottleneck, while voxel-pair filtering completes quickly. Here, the advantage comes from our best-first search strategy, which prioritizes candidates with smaller distances and prunes unpromising pairs earlier, whereas TDBase adopts a less efficient depth-first traversal. Refinement Stage. Figure 16 shows the refinement-stage speedup of 3DPipe over TDBase for both within-𝜏 and 𝑘-NN queries across all datasets, where the x-axis varies 𝑘 or 𝜏. Overall, 3DPipe achieves up to 4× speedup on NV/NN (Figure 16a) and up to 10× on TI/TT (Figure 16b). Notably, although TDBase already leverages GPU acceleration, its performance remains suboptimal due to inefficient hardware utilization. In contrast, 3DPipe consistently delivers higher performance by fully exploiting GPU capabilities. In particular, workload flattening

Conference acronym ’XX, June 03–05, 2018, Woodstock, NY

k-NN TI Within-τ TI k-NN TT Within-τ TT

k-NN NV Within-τ NV k-NN NN Within-τ NN

k or τ

50 0 Within-400 10-NN Within-200 10-NN Within-20 10-NN Within-10 50-NN

(b) Refinement on TI and TT

Figure 16: Runtime Speedup of Refinement Stage Unified Memo y

150 50 30 40 100 20 200 30 50 A100 Memory Capacity 20 A100 Memory Capacity 50 10 10 100 0 0 50 100150200250300350400 0 1 2 3 4 5 6 7 8 9 10 0 k (a) Within-𝜏 𝜏) (b) 𝑘-NN NV (varying 𝑘) 0 1 2 3NV4(varying 5 6 7 8 9 10 K 1016 Unified Memor Unified Memory

0 Within-400 10-NN Within-200 10-NN Within-20 10-NN Within-10 50-NN NV

GPU

109 Algorithm 1

Algorithm 2

(d) 10-NN NV

Figure 17: Unified Memory v.s. Chunked Streaming

10 1

300 200

0.1 1 2 3 4 5 6100 7 8 9 10 k

(a) 𝑘-NN NV (varying 𝑘)

and shared-memory aggregation significantly improve parallel efficiency, while CPU-GPU pipelining overlaps data preparation with computation to further enhance utilization.

4.4

Effectiveness of Chunked Streaming

Recall that our chunked streaming design (Algorithms 3 and 5) enables processing voxel pairs beyond GPU memory capacity. We compare it with unified memory (cudaMallocManaged), which relies on on-demand page migration so can exceed the 40GB A100 memory capacity. Figure 17 reports the performance on NV for both within-𝜏 and 𝑘-NN queries, where the curves show runtime (left y-axis), and the bars indicate memory usage (right y-axis). In Figure 17a, when 𝜏 ≤ 250, the working set fits within A100 memory, and our approach already outperforms unified memory with a modest margin. Once 𝜏 ≥ 300, memory demand exceeds capacity, and unified memory exhibits a sharp slowdown, while our approach remains stable, leading to a much larger performance gap. This effect is more pronounced in Figure 17b, since memory demand exceeds capacity even at 𝑘 = 1. To further understand this gap, we profile two representative queries (within-400 and 10-NN) using NVIDIA Nsight Compute. Figures 17c and 17d report GPU memory stall cycles, which indicate stalled warps due to memory access delays. Unified memory incurs 146× and 184× more stall cycles than our approach for Algorithms 1 and 2, respectively. These results show that unified memory suffers from significant overhead once page migration is triggered, whereas

Figure 18: CPU-GPU Pipelining v.s. No Pipelining

Runtime (sec)

Total GPU Memory Stall Cycles

Algorithm 2

CPU–GPU Pipelining

100

1012

109 Algorithm 1

(b) GPU Memory Consumption of Refinement Stage

Ch nked Streaming Approach

1014

1012

Space (GB)

Runtime (sec) Runtime (sec) Total GPU Memor Stall C cles

Chunked Streaming Approach

1014

No Pipeline

Runtime (sec)

100

300

(a) Running Time of Refinement Stage

GPU Memor (GB)

Chunked St eaming App oach

CPU–GPU Pipelining

100

1 2 3 4 5 6 7 8 9 10 k or τ (a) Refinement on NV and NN

Data Preparation (No pipe) Compute (No pipe)

CPU 100 Runtime (sec)

Speedup

Runtime ( )

Lyuheng Yuan, Da Yan, Akhlaque Ahmad and Fusheng Wang

10 1 1 2 3 4 5 6 7 8 9 10 k

(b) 𝑘-NN NN (varying 𝑘)

1 2 3 4 Pruning 5 6 7 by8CPU 9 10 Figure 19: 𝑘-NN0Object-Pair v.s. GPU K

chunked streaming avoids this issue, achieving up to 5.4× and 8.9× speedup for within-400 and 10-NN queries, respectively.

4.5

Effectiveness of CPU-GPU Pipelining

Recall that refinement involves both CPU data preparation and GPU computation. We compare our CPU-GPU pipelining (Algorithm 5) with a non-pipelined baseline, where CPU data preparation and GPU computation are executed sequentially. Figure 18 shows the impact on runtime and GPU memory usage. As shown in Figure 18a, pipelining reduces the refinement time by up to 69% by overlapping CPU and GPU execution. It also significantly reduces memory usage: since at most two chunks reside on the GPU simultaneously, the peak memory consumption is only 1.7 GB, compared to 16.5 GB for the baseline (Figure 18b), while eliminating GPU idle time during data preparation.

4.6

Evaluation of Object-Pair Pruning on GPU

Recall that for a 𝑘-NN query, object-pair pruning is repeatedly invoked after filtering and each refinement round. TDBase performs this step on CPU, whereas we design a GPU-based kernel (Algorithm 6) to fully exploit parallelism. Figure 19 compares the total runtime of object-pair pruning (summed over all its invocations). Overall, the GPU-based approach achieves nearly two orders of magnitude speedup over the CPU

3DPipe: GPU-Accelerated 3D Spatial Join

Conference acronym ’XX, June 03–05, 2018, Woodstock, NY

Pipeline

No Pipeline

Figure 21: Screenshots from Nvidia Nsight Systems. Top: Pipelining with CUDA Stream. Bottom: No Pipelining.

400

End-to-End Time (sec)

Runtime (sec)

w/ Pipeline w/o Pipeline

10-NN NV Within-400 NV 10-NN NN Within-200 NN

200

9 10

Global Memory Shared Memory

Runtime (sec)

300 200

Total GPU Memory Stall Cycles

Figure 20: Time of Algorithm 3 With v.s. Without Pipeline 1014

0 Within-400 NV 10-NN NV Within-200 NN 10-NN NN (a) Runtime of Refinement

109 Within-400 NV 10-NN NV Within-200 NN 10-NN NN (b) Memory Stall Cycles

Figure 22: Aggregation in GPU: Global v.s. Shared Memory baseline with OpenMP. On NV (Figure 19a), it is 23×–99× faster, where CPU execution is dominated by expensive filtering rounds (up to 42 s), while the GPU completes each round in 0.04 s. On NN (Figure 19b), the speedup remains 9×–42×. In the final refinement round (LoD = 100%), where CPU-based refinement dominates (up to 49 s), the GPU still completes selection in about one second. These results show that GPU-based object-pair pruning is critical to performance and scales well with the number of candidates.

4.7

Other Ablation Studies

We present additional ablation studies on secondary optimizations, complementing the major components evaluated above. CUDA Stream Pipeline. Our chunked streaming handles data beyond GPU memory capacity, but GPU may still remain idle during device-to-host data transfer. We introduce a pipelined implementation (Algorithm 3) to overlap computation and data transfer. Figure 20 compares chunked streaming with and without pipelining on 𝑘-NN queries over NV. We observe a consistent ∼10% speedup, as pipelining utilizes idle GPU cycles by initiating computation for the next chunk while transferring results. To verify this overlap, we profile execution using NVIDIA Nsight Systems. As Figure 21 shows, kernel execution overlaps with memory transfer when pipelining is enabled, whereas no overlap is observed otherwise. This confirms that asynchronous execution effectively hides data transfer latency and improves GPU utilization.

8 |V| (x103)

Figure 23: End-to-End Time with Increasing Data Scale

Global Memory Shared Memory

1012

100

Aggregation with Shared Memory. Recall that block-wise aggregation (Figure 6) computes the aggregation across a thread block using shared memory, which is used in both filtering and refinement of 3DPipe. In particular, our refinement kernel, Algorithm 4, uses it in Line 13. TDBase instead performs this minimum-aggregation in global memory using atomicMin. Figure 22 quantifies the impact of this shared-memory aggregation on refinement performance. As shown in Figure 22a, sharedmemory aggregation achieves 6×–7× speedup over global-memory aggregation. Meanwhile, Figure 22b reports GPU memory stall cycles (caused by cache/global memory accesses). Shared memory reduces stall cycles by 32×–108×, indicating significantly lower memory latency and improved warp scheduling efficiency.

4.8

Scalability Analysis

To evaluate scalability, we generate vessel and nuclei datasets with increasing sizes from 1K to 16K, doubling at each step. We use four representative queries with the large workloads. Figure 23 shows that our approach scales nearly linearly with data size, while maintaining high efficiency across all settings. This demonstrates its ability to handle large-scale datasets effectively.

Related Work

Spatial Join Processing. Spatial join has been studied extensively for 2D spatial data. Early work mainly focused on disk-based methods, such as hash-based [27] and tile-based spatial join [9, 35], to cope with limited main memory. As memory capacity increased, in-memory methods [33] and later parallel and distributed systems [8, 10, 12, 16, 47, 48, 52] were developed to improve performance on large-scale datasets. In addition to CPU-based acceleration, modern hardware such as GPU and FPGA has also been explored for faster spatial join processing [22, 50]. Mesh Simplification and Multi-LoD Representations. To reduce the cost of geometric computation in 3D, prior work has

Conference acronym ’XX, June 03–05, 2018, Woodstock, NY

explored progressive mesh simplification [24, 43], which constructs multiple levels of detail (LoDs) for the same polyhedron. Such multi-LoD representations are useful for spatial query processing because they enable coarse-to-fine evaluation: many object pairs can be pruned or resolved at low resolutions, avoiding unnecessary computation on the original high-resolution geometry. 3D Spatial Data Management and Query Processing. Compared with the rich literature on 2D spatial join, 3D spatial join over polyhedral objects remains much less explored. Among existing systems, iSPEED [25] adopts a MapReduce-based framework for distributed 3D spatial queries, with optimizations to reduce disk I/O. Real et al. [36] accelerate basic 3D spatial operations on top of PostGIS [6]. However, general-purpose spatial databases such as PostGIS lack efficient support for large-scale 3D spatial join, especially when polyhedra contain many facets and complex geometry, which motivates specialized solutions. More closely related to our work, 3DPro [39] introduces progressive refinement over multiple LoDs and leverages GPU to accelerate the refinement stage. TDBase [40] further improves this framework by incorporating facet-level Hausdorff and proxy Hausdorff bounds, enabling tighter distance estimation and more effective pruning at low LoDs.

Conclusions and Future Work

We presented 3DPipe, a pipelined GPU framework for scalable 3D spatial join over polyhedral objects. Our design exploited GPU parallelism across both filtering and refinement, and introduced chunked streaming and CPU-GPU pipelining to overcome memory and utilization bottlenecks. Together with GPU-oriented optimizations such as shared-memory aggregation and parallel top-𝑘 filtering, 3DPipe achieved efficient end-to-end execution that consistently outperformed the state-of-the-art solution TDBase, achieving up to 9.0× speedup while maintaining near-linear scalability. Currently, all facet data are loaded in host memory, extracted and transferred to the GPU on demand during refinement. While effective, this design incurs non-trivial host memory overhead for large datasets. In future work, we plan to extend 3DPipe to an out-of-core setting by storing facet data on SSD and adopting a index nested-loop join strategy. Given that facet-level distance computation is the dominant cost and is highly GPU-parallelizable, we expect that pipelined execution can overlap SSD I/O with GPU computation and effectively hide I/O latency.

References [1] [n. d.]. Cesium. https://cesium.com. [2] [n. d.]. ERSI 3D GIS. https://www.esri.com/en-us/capabilities/3d-gis. [3] [n. d.]. ModelNet40 - Princeton 3D Object Dataset. https://www.kaggle.com/dat asets/balraj98/modelnet40-princeton-3d-object-dataset. [4] [n. d.]. OFF (Object File Format). https://en.wikipedia.org/wiki/OFF_(file_format ). [5] [n. d.]. Polyhedron Model in PolarDB. https://www.alibabacloud.com/help/en/p olardb/polardb-for-postgresql/models-pg. [6] [n. d.]. PostGIS. https://postgis.net/. [7] [n. d.]. TDBase Source Code. https://github.com/tengdj/tdbase. [8] Ablimit Aji, Fusheng Wang, Hoang Vo, Rubao Lee, Qiaoling Liu, Xiaodong Zhang, and Joel Saltz. 2013. Hadoop-GIS: A high performance spatial data warehousing system over MapReduce. In Proceedings of the VLDB endowment international conference on very large data bases, Vol. 6. p1009. [9] Lars Arge, Octavian Procopiuc, Sridhar Ramaswamy, Torsten Suel, and Jeffrey Scott Vitter. 1998. Scalable sweeping-based spatial join. In VLDB, Vol. 98. 570–581.

Lyuheng Yuan, Da Yan, Akhlaque Ahmad and Fusheng Wang

[10] Furqan Baig, Hoang Vo, Tahsin Kurc, Joel Saltz, and Fusheng Wang. 2017. Sparkgis: Resource aware efficient in-memory spatial query processing. In Proceedings of the 25th ACM SIGSPATIAL international conference on advances in geographic information systems. 1–10. [11] Fausto Bernardini, Joshua Mittleman, Holly E. Rushmeier, Cláudio T. Silva, and Gabriel Taubin. 1999. The Ball-Pivoting Algorithm for Surface Reconstruction. IEEE Trans. Vis. Comput. Graph. 5, 4 (1999), 349–359. [12] Thomas Brinkhoff, H-P Kriegel, and Bernhard Seeger. 1996. Parallel processing of spatial joins using R-trees. In Proceedings of the Twelfth International Conference on Data Engineering. IEEE, 258–265. [13] Holger Caesar, Varun Bankiti, Alex H. Lang, Sourabh Vora, Venice Erin Liong, Qiang Xu, Anush Krishnan, Yu Pan, Giancarlo Baldan, and Oscar Beijbom. 2020. nuScenes: A Multimodal Dataset for Autonomous Driving. In CVPR. Computer Vision Foundation / IEEE, 11618–11628. [14] Xingyu Chen, Fu-Jen Chu, Pierre Gleize, Kevin J. Liang, Alexander Sax, Hao Tang, Weiyao Wang, Michelle Guo, Thibaut Hardin, Xiang Li, Aohan Lin, Jiawei Liu, Ziqi Ma, Anushka Sagar, Bowen Song, Xiaodong Wang, Jianing Yang, Bowen Zhang, Piotr Dollár, Georgia Gkioxari, Matt Feiszli, and Jitendra Malik. 2025. SAM 3D: 3Dfy Anything in Images. arXiv preprint arXiv:2511.16624 (2025). doi:10 .48550/arXiv.2511.16624 [15] HuBMAP Consortium. 2019. The human body at cellular resolution: the NIH Human Biomolecular Atlas Program. Nature 574, 7777 (2019), 187–192. [16] Ahmed Eldawy and Mohamed F Mokbel. 2015. Spatialhadoop: A mapreduce framework for spatial data. In 2015 IEEE 31st international conference on Data Engineering. IEEE, 1352–1363. [17] Navid Farahani, Alex Braun, Dylan Jutt, Todd Huffman, Nick Reder, Zheng Liu, Yukako Yagi, and Liron Pantanowitz. 2017. Three-dimensional imaging and scanning: current and future applications for pathology. Journal of pathology informatics 8, 1 (2017), 36. [18] Andreas Geiger, Philip Lenz, and Raquel Urtasun. 2012. Are we ready for autonomous driving? The KITTI vision benchmark suite. In CVPR. IEEE Computer Society, 3354–3361. [19] S Gnanakaran, Hugh Nymeyer, John Portman, Kevin Y Sanbonmatsu, and Angel E Garcıa. 2003. Peptide folding simulations. Current opinion in structural biology 13, 2 (2003), 168–174. [20] Ziyu Guo, Renrui Zhang, Xiangyang Zhu, Chengzhuo Tong, Peng Gao, Chunyuan Li, and Pheng-Ann Heng. 2024. SAM2Point: Segment Any 3D as Videos in Zeroshot and Promptable Manners. CoRR abs/2408.16768 (2024). [21] W. Daniel Hillis and Guy L. Steele Jr. 1986. Data Parallel Algorithms. Commun. ACM 29, 12 (1986), 1170–1183. doi:10.1145/7902.7903 [22] Wenqi Jiang, Oleh-Yevhen Khavrona, Martin Parvanov, and Gustavo Alonso. 2025. Swiftspatial: Spatial joins on modern hardware. Proceedings of the ACM on Management of Data 3, 3 (2025), 1–27. [23] Michael M. Kazhdan, Matthew Bolitho, and Hugues Hoppe. 2006. Poisson surface reconstruction. In Proceedings of the Fourth Eurographics Symposium on Geometry Processing (ACM International Conference Proceeding Series), Alla Sheffer and Konrad Polthier (Eds.). Eurographics Association, 61–70. [24] Andrei Khodakovsky, Peter Schröder, and Wim Sweldens. 2000. Progressive geometry compression. In Proceedings of the 27th annual conference on Computer graphics and interactive techniques. 271–278. [25] Yanhui Liang, Hoang Vo, Jun Kong, and Fusheng Wang. 2017. iSPEED: an Efficient In-Memory Based Spatial Query System for Large-Scale 3D Data with Complex Structures. In Proceedings of the 25th ACM SIGSPATIAL International Conference on Advances in Geographic Information Systems, GIS 2017, Redondo Beach, CA, USA, November 7-10, 2017. ACM, 17:1–17:10. [26] Yanhui Liang, Fusheng Wang, Pengyue Zhang, Joel H Saltz, Daniel J Brat, and Jun Kong. 2017. Development of a framework for large scale three-dimensional pathology and biomarker imaging and spatial analytics. AMIA Summits on Translational Science Proceedings 2017 (2017), 75. [27] Ming-Ling Lo and Chinya V Ravishankar. 1996. Spatial hash-joins. In Proceedings of the 1996 ACM SIGMOD international conference on Management of data. 247– 258. [28] William E. Lorensen and Harvey E. Cline. 1987. Marching cubes: A high resolution 3D surface construction algorithm. In SIGGRAPH, Maureen C. Stone (Ed.). ACM, 163–169. [29] Adrien Maglo, Clement Courbet, Pierre Alliez, and Céline Hudelot. 2012. Progressive compression of manifold polygon meshes. Comput. Graph. 36, 5 (2012), 349–359. [30] Adrien Maglo, Guillaume Lavoué, Florent Dupont, and Céline Hudelot. 2015. 3D Mesh Compression: Survey, Comparisons, and Emerging Trends. ACM Comput. Surv. 47, 3 (2015), 44:1–44:41. [31] Ben Mildenhall, Pratul P. Srinivasan, Matthew Tancik, Jonathan T. Barron, Ravi Ramamoorthi, and Ren Ng. 2020. NeRF: Representing Scenes as Neural Radiance Fields for View Synthesis. In ECCV (Lecture Notes in Computer Science), Andrea Vedaldi, Horst Bischof, Thomas Brox, and Jan-Michael Frahm (Eds.). Springer, 405–421. [32] Tomas Möller. 1997. A Fast Triangle-Triangle Intersection Test. J. Graphics, GPU, & Game Tools 2, 2 (1997), 25–30.

3DPipe: GPU-Accelerated 3D Spatial Join

[33] Sadegh Nobari, Farhan Tauheed, Thomas Heinis, Panagiotis Karras, Stéphane Bressan, and Anastasia Ailamaki. 2013. TOUCH: in-memory spatial join by hierarchical data-oriented partitioning. In Proceedings of the ACM SIGMOD International Conference on Management of Data, SIGMOD 2013, New York, NY, USA, June 22-27, 2013. ACM, 701–712. [34] Jeong Joon Park, Peter R. Florence, Julian Straub, Richard A. Newcombe, and Steven Lovegrove. 2019. DeepSDF: Learning Continuous Signed Distance Functions for Shape Representation. In CVPR. Computer Vision Foundation / IEEE, 165–174. [35] Jignesh M. Patel and David J. DeWitt. 1996. Partition Based Spatial-Merge Join. In Proceedings of the 1996 ACM SIGMOD International Conference on Management of Data, Montreal, Quebec, Canada, June 4-6, 1996. ACM Press, 259–270. [36] Lucas C Villa Real, Bruno Silva, Dikran S Meliksetian, and Kaique Sacchi. 2019. Large-scale 3D geospatial processing made possible. In Proceedings of the 27th ACM SIGSPATIAL International Conference on Advances in Geographic Information Systems. 199–208. [37] Nick Roussopoulos, Stephen Kelley, and Frédéic Vincent. 1995. Nearest Neighbor Queries. In SIGMOD, Michael J. Carey and Donovan A. Schneider (Eds.). ACM Press, 71–79. [38] Dejun Teng, Furqan Baig, Zhaohui Peng, Jun Kong, and Fusheng Wang. 2024. Efficient spatial queries over complex polygons with hybrid representations. GeoInformatica 28, 3 (2024), 459–497. [39] Dejun Teng, Furqan Baig, Hoang Vo, Yanhui Liang, Jun Kong, and Fusheng Wang. 2022. 3DPro: Querying Complex Three-Dimensional Data with Progressive Compression and Refinement. In Proceedings of the 25th International Conference on Extending Database Technology, EDBT 2022, Edinburgh, UK, March 29 - April 1, 2022. 2:104–2:117. [40] Dejun Teng, Zhaochuan Li, Zhaohui Peng, Shuai Ma, and Fusheng Wang. 2025. Efficient and Accurate Spatial Queries Using Lossy Compressed 3D Geometry Data. IEEE Trans. Knowl. Data Eng. 37, 5 (2025), 2472–2487. [41] Dejun Teng, Yanhui Liang, Hoang Vo, Jun Kong, and Fusheng Wang. 2022. Efficient 3D Spatial Queries for Complex Objects. ACM Trans. Spatial Algorithms Syst. 8, 2 (2022), 1–26. [42] Haithem Turki, Deva Ramanan, and Mahadev Satyanarayanan. 2022. Mega-NeRF: Scalable Construction of Large-Scale NeRFs for Virtual Fly-Throughs. In CVPR. IEEE, 12912–12921.

Conference acronym ’XX, June 03–05, 2018, Woodstock, NY

[43] Sébastien Valette, Raphaëlle Chaine, and Rémy Prost. 2009. Progressive lossless mesh compression via incremental parametric refinement. In Computer Graphics Forum, Vol. 28. 1301–1310. [44] Fusheng Wang, Jun Kong, Lee Cooper, Tony Pan, Tahsin Kurc, Wenjin Chen, Ashish Sharma, Cristobal Niedermayr, Tae W Oh, Daniel Brat, et al. 2011. A data model and database for high-resolution pathology analytical image informatics. Journal of pathology informatics 2, 1 (2011), 32. [45] World Labs. [n. d.]. Marble Labs: Blueprints for Building with World Models. https://www.worldlabs.ai/labs. [46] Zhirong Wu, Shuran Song, Aditya Khosla, Fisher Yu, Linguang Zhang, Xiaoou Tang, and Jianxiong Xiao. 2015. 3d shapenets: A deep representation for volumetric shapes. In Proceedings of the IEEE conference on computer vision and pattern recognition. 1912–1920. [47] Jia Yu, Jinxuan Wu, and Mohamed Sarwat. 2015. GeoSpark: a cluster computing framework for processing large-scale spatial data. In Proceedings of the 23rd SIGSPATIAL International Conference on Advances in Geographic Information Systems, Bellevue, WA, USA, November 3-6, 2015. ACM, 70:1–70:4. [48] Lyuheng Yuan, Da Yan, Akhlaque Ahmad, Jiao Han, Saugat Adhikari, and Yang Zhou. 2025. Out-of-Core Parallel Spatial Join Outperforming In-Memory Systems: A BFS-DFS Hybrid Approach. In Proceedings of the 34th International Symposium on High-Performance Parallel and Distributed Computing, HPDC 2025, Notre Dame, IN, USA, July 20-23, 2025. ACM, 23:1–23:14. [49] Andi Zang, Shiyu Luo, Xin Chen, and Goce Trajcevski. 2019. Real-Time Applications Using High Resolution 3D Objects in High Definition Maps (Systems Paper). In SIGSPATIAL. ACM, 229–238. [50] Jianting Zhang, Simin You, and Le Gruenwald. 2017. Parallel selectivity estimation for optimizing multidimensional spatial join processing on gpus. In 2017 IEEE 33rd International Conference on Data Engineering (ICDE). IEEE, 1591–1598. [51] Rui Zhang, Jianzhong Qi, Dan Lin, Wei Wang, and Raymond Chi-Wing Wong. 2012. A highly optimized algorithm for continuous intersection join queries over moving objects. VLDB J. 21, 4 (2012), 561–586. [52] Xiaofang Zhou, David J Abel, and David Truffet. 1998. Data partitioning for parallel spatial join processing. Geoinformatica 2, 2 (1998), 175–204. [53] Yuchen Zhou, Jiayuan Gu, Tung Yen Chiang, Fanbo Xiang, and Hao Su. 2025. Point-SAM: Promptable 3D Segmentation Model for Point Clouds. In ICLR. OpenReview.net.

3DPipe: A Pipelined GPU Framework for Scalable Generalized Spatial Join over Polyhedral Objects

Related documents