Spatial join queries are pivotal in Geographic Information Science (GIS) and geo-databases. A spatial join query aims to identify pairs of geometries from multiple datasets that satisfy a given predicate, such as "intersect" or "contain. " [7,32,45]. For example, it can determine the locations of bridges by querying intersections between a road network and a hydrographic map. However, the increasing volume of GIS data collected from various sources, including satellites and mobile devices, poses challenges for timely data analytics [4,10,37,51]. This requirement is particularly acute in domains like interactive urban planning, emergency management, and risk assessment [58,72,79].
Various techniques have been explored to address the challenges of large-scale data and low-latency processing. The plane-sweep algorithm is favored for its low complexity [12,47,66] but suffers from limited parallel processing capabilities [5,6,41]. Rasterizationbased methods offer parallelism but require intensive preprocessing [16,65,79]. Learned spatial indexes, though efficient in query performance and storage, struggle with high training costs and accuracy concerns [18,29,52,53,59,82]. Consequently, GIS researchers and practitioners often resort to conventional indexing techniques, such as grids and trees, which offer numerous advantages [4, 6, 10, 17, 62-64, 71, 73-77, 80].
However, the conventional spatial indexes mentioned above have performance limitations from achieving real-time processing [61]. While tree-based methods offer a unique advantage in logarithmic time complexity, they suffer from inherent shortcomings, including branches, irregular memory accesses, and memory space overhead for storing pointers. These drawbacks notably undermine the efficiency of spatial join operations, particularly on GPUs. In contrast, grids convert spatial queries into hardware-friendly scan operations within a limited search space, thereby maximizing memory access efficiency. However, skewed data distributions may result in imbalanced workloads that hinder the effectiveness of parallel processing [43,44,69].
Over the years, algorithmic improvements in conventional indexing methods seriously mismatch with the execution models of general-purpose computer architecture, including CPUs and GPUs, leading to persistent performance bottlenecks and associated tradeoffs in designing databases and analytical systems [13,31,78,81]. In fact, our evaluation shows that none of the existing spatial join methods satisfy the real-time processing requirement ( §5.7), defined as the completion time within a second. Recently, the GPU advancement of dedicated ray-tracing units, such as NVIDIA's RT Cores, opens a new opportunity to accelerate spatial join queries through hardware-accelerated Bounding Volume Hierarchy (BVH) tree traversal. Moreover, the processing capability of RT Cores is doubled per generation, making RT-based methods automatically benefit from the advancing of this new hardware [48,49]. Utilizing the technology advancement, we have designed and implemented a new and effective acceleration engine for spatial joins with RT Cores, which is called RayJoin, to address the limitations of all existing approaches.
Leveraging RT Cores for spatial join queries presents three nontrivial challenges for us to address. First, to harness the dedicated hardware effectively, it is crucial to devise an efficient method that formulates spatial join queries as ray tracing problems. In RayJoin, we showcase two critical spatial join queries, line segment intersection (LSI) and point-in-polygon (PIP) test. We transform these two queries into RT-friendly algorithms using the 𝐴𝑛𝑦𝐻𝑖𝑡 and 𝐶𝑙𝑜𝑠𝑒𝑠𝑡𝐻𝑖𝑡 shaders supported by RT Cores. Specifically, we have leveraged these queries to develop a polygon overlay application. Second, the precision in RT Cores is limited to single-precision floating-point (FP32) for high-performance rendering. However, this limitation poses a significant challenge in meeting the requirements of GIS applications. To address this, we have proposed a conservative representation that stores polygons in low precision while still delivering exact query results. In computer arithmetic, a conservative method ensures the precision of computation results, even at a cost of performance. RayJoin achieves both high precision and high performance. Finally, the BVH construction cost, in both time and space, is another bottleneck to be resolved to achieve high performance. We have noticed that the construction time is often a magnitude higher than the query time. In response, we have introduced a grouping technique that greatly reduces the construction cost by reducing the number of primitives used in the BVH. This effectively reduces construction costs while keeping high query performance.
Our experiments show that RayJoin has achieved 3.0x to 28.3x speedups over any existing polygon overlay solutions. It only takes about 460ms on the full-fledged polygon overlay analysis to join 1.6M polygons with 303K polygons, meeting real-time spatial join query requirements. In this paper, we have made a strong case for using RT Cores to accelerate various spatial queries in geodatabases, and our contributions are as follows:
• We have developed an effective approach that utilizes RT Cores to accelerate LSI and PIP queries, substantially outperforming conventional spatial join techniques on both CPUs and GPUs. To the best of our knowledge, RayJoin is the fastest spatial join solution ( §3.1). • Within the existing GPUs that only support FP32 ray-tracing, we have developed a method to effectively achieve doubleprecision floating-point (FP64) while maintaining performance levels similar to the underlying hardware ( §3.2). • We have proposed an algorithm that greatly reduces the BVH construction time and memory cost while keeping high query performance ( §3.3). • We have comprehensively compared and analyzed RayJoin with many existing spatial join solutions by a consistent benchmark with synthetic and real-world datasets ( §5), making a strong case for RT Cores to attain fast and precise spatial join. • RayJoin is an open-source software project for widespread public usage.
1 1 Source code: https://github.com/pwrliang/RayJoin 2 BACKGROUND 2.1 Spatial Join Queries and Spatial Indexes Spatial join finds pairs of spatial objects (𝑟, 𝑠) from two datasets 𝑅 and 𝑆 that satisfy a spatial predicate 𝜃 , where 𝑟 ∈ 𝑅, 𝑠 ∈ 𝑆. Formally, 𝑅 ⊲⊳ 𝜃 𝑆 {(𝑟, 𝑠) | 𝑟 ∈ 𝑅, 𝑠 ∈ 𝑆, 𝜃 (𝑟, 𝑠) = 𝑡𝑟𝑢𝑒}
Spatial objects come in several forms, including points, line segments, polylines, and polygons. Predicates, such as intersect, contain, overlap, etc., can be combined to perform complex spatial queries. Figure 1 provides an example of joining two sets of polygons 𝑅 and 𝑆, which is also called polygon overlay analysisfoot_0 . The polygon overlay is achieved by combining the results of LSI and PIP queries. The LSI query generates points at the intersections where the boundaries of datasets 𝑅 and 𝑆 meet. The PIP query determines the polygon 𝑟 ∈ 𝑅 where a query point 𝑠 ∈ 𝑆 lies and vice versa. The output from the above two steps is combined to obtain the resulting polygon (on the right of Figure 1).
Performing spatial join queries on a large volume of complex geometries is challenging due to the high computational cost [3]. Exhaustive search compares all pairs of geometries, resulting in 𝑂 (|𝑅| • |𝑆 |) time complexity, which is impractical for large datasets [40]. Thus, spatial indexes are necessary to improve query performance by reducing search space. This section introduces two typical tree-based and grid-based spatial indexes, namely BVH and uniform grid. To clarify, in graphics, BVH is a general term to describe a hierarchical data structure that partitions geometries with bounding volume. For example, the R-tree can be categorized into BVH because it uses Minimum Bounding Rectangle (MBR) as the bounding volume [19].
BVH. The BVHs organize geometries in tree-like data structures. The internal node in the BVH is associated with an axis-aligned bounding box (AABB) 3 , which tightly encloses the union of the children's AABBs. Leaf nodes hold the primitives 4 , which can be points, lines, triangles, or any shapes. AABB approximates the primitive boundary, making the intersection test more efficient. Building a BVH can be efficiently converted into a GPU-friendly
Logical pointers sorting problem. For example, Linear Bounding Volume Hierarchy (LBVH) organizes geometries in the order of their Morton Code, balancing the construction time and index quality [26].
Figure 2 (a) depicts a query 𝑄 𝐸 being performed on a scene composed of four triangles. Figure 2 (b) illustrates how to accelerate the query with the BVH. The query aims to identify all the triangles intersecting the query. This process begins by performing an intersection test between 𝑄 𝐸 and the AABB 𝐵 1 of the root node. If an intersection is detected, further tests are conducted between 𝑄 𝐸 and the nodes 𝑁 2 and 𝑁 3 . The search proceeds until we narrow it down to the leaf node 𝑁 2 . Then, we only need to execute the intersection test against the triangle bounded by 𝐵 2 . The right sub-tree rooted at 𝑁 3 can then be pruned as it does not intersect with the query.
Tree-based indexes offer logarithmic time complexity of searching, which is a key advantage. However, they can be memory inefficient. As shown in Figure 2 (b), tree nodes are stored in the "𝑁𝑜𝑑𝑒𝑠" array. Logical pointers connect parent nodes to their child nodes. Each query starts by visiting the root node and following the logical pointers to its descendants, resulting in a Depth-First Search (DFS). The pointer-chasing during DFS leads to irregular memory access patterns and branches, preventing coalesced memory access optimization and resulting in cache line inefficiencies. As a result, numerous memory transactions are generated during the traversal. Determining the traversal path also results in branch divergence on the GPU, causing very few threads to be active in each warp.
Uniform Grid. The uniform grid is a spatial index that partitions the space into equal-sized cells. The width of the grid is called the resolution. Figure 2 (c) represents a 3x3 uniform grid built for the scene. The triangles are split into line segments for efficient retrieval. The grid data structure consists of "Cells" and "Line Segs" array. The "Cells" array maintains the beginning offsets that point to the positions in the "Line Segs" array, where the line segments within the cell are organized in a contiguous memory region.
Accelerating the query 𝑄 𝐸 can be achieved through the uniform grid, as demonstrated in Figure 2 (c). The cells that overlap with the query are collected (shaded with gray). Then, we scan line segments in the cells to determine whether they intersect with the query. With the grid, the line segments spatially distant from the query are pruned. The finer the grid, the more line segments are pruned, leading to improved query efficiency, but the construction cost is also increased due to duplicated line segments. Thus, the resolution of the grid must be carefully tuned.
The uniform grid offers several advantages, including GPUfriendly scan operations and a low construction cost. If geometries are distributed uniformly, the grid-based indexes should be ideal candidates for GPUs. Unfortunately, this is not true in real-world applications. For example, geometries representing buildings are dense in metropolitan cities but usually sparse in deserts, which forms an uneven grid, limiting the parallelism of GPUs and leading to load imbalance and low occupancy.
Ray tracing is a technique for rendering photorealistic images. The fundamental principle of ray tracing is identifying the primitives hit by rays, which is time-consuming, typically requiring an Accelerating Structure (AS) to reduce search space. Currently, BVH is the most widely utilized AS in ray tracing frameworks, as discussed in §2.1. Despite their efficiency, tree-based indexes are not optimal for use on GPUs. As a result, researchers have proposed various hardware-based solutions to accelerate BVH traversal. For example, NVIDIA RTX series GPUs provide dedicated ray tracing units to accelerate both BVH traversal and ray-primitive intersection test [8]. Several ray-tracing interfaces are available for the ray-tracing hardware, including NVIDIA OptiX, Vulkan, and DirectX Raytracing. This paper focuses on OptiX as it is interoperable with CUDA.
OptiX supports the "Multiple Instruction, Multiple Data" (MIMD) subset of CUDA [50]. Each thread casts a ray, and the execution of the ray could be rescheduled at runtime to a different lane, warp, and even Streaming Multiprocessor (SM) for execution efficiency. Thus, shared memory and warp/block synchronization intrinsics are not available. A ray is allowed to carry per-ray data via registers to communicate to shaders. Shaders are user-defined CUDA programs to cast rays and handle ray-primitive intersection events.
Before casting rays, we have to build the BVH from an array of primitives. The BVH structure is hardware-specific and handled by the video driver. Thus, users have no control over the BVH construction process but only supply primitives. Constructing a BVH is expensive. The construction time is linearly correlated to the number of primitives [83]. Upon completion of the BVH construction, a traversal handle is returned, which is used to access the BVH. The BVH traversal occurs on the RT Cores, which share the same device memory with the SMs.
Figure 3 illustrates the execution flow of OptiX (shapes shaded in grey are non-programmable procedures). OptiX does not allow explicit BVH traversal and instead begins by generating a set of rays. The entry point of a ray tracing program is the RayGen (RG) shader, where the traversal handle, origins, and directions of the rays are provided to generate rays and identify hits. The BVH is
SMs RT Cores TL 1 TL 2 RayGen (RG) Traversal Handle Ray primitive intersect? BVH Traversal (TL) IsIntersection (IS) AnyHit (AH) Found a hit? ClosestHit (CH) Miss (MS) CH IS RG Yes No Traversal completes Yes No Leaf node TL 3 AH Figure 3: Top: Switching between the RT Cores and SMs; Bottom: Execution flow of OptiX. autonomously traversed, guided by the rays. The traversal stage is non-programmatic and transparent to the user. Figure 3 top shows the execution flow of the query on Figure 2 (b), where 𝑇 𝐿 𝑖 corresponds to the traversal of BVH node 𝑁 𝑖 . If leaf nodes are reached during the traversal, the IsIntersection (IS) shader is invoked to determine if a ray intersects with a primitive. When the IS shader reports an intersection, the AnyHit (AH) shader is called, allowing the user to execute rendering tasks. At the end of the traversal, the ClosestHit (CH) shader is called, which returns the closest primitive hit by the ray. The Miss (MS) shader is invoked if the ray does not hit any primitive. In OptiX, the coordinates of primitives and ray parameters are defined in the FP32 for performance reasons.
This limitation makes implementing high-precision applications challenging. Thus, this issue must be addressed to accommodate broad applications.
Ray Definition. A ray R is a semi-infinite line characterized by its origin 𝑂, a direction vector ì 𝑑, and a parameter 𝑡. The parameter 𝑡 allows the ray to extend indefinitely in the direction of ì 𝑑 starting from 𝑂. Formally, it is defined in Equation (1):
When casting a ray, users are asked to provide a range parameter [𝑡 𝑚𝑖𝑛 , 𝑡 𝑚𝑎𝑥 ] to filter the intersections on a segment of the ray (Figure 4 (a)). If the ray hits a primitive (such as an AABB), a hit event is reported only if 𝑡 𝑚𝑖𝑛 ≤ 𝑡 ℎ𝑖𝑡 ≤ 𝑡 𝑚𝑎𝑥 , where 𝑂 + 𝑡 ℎ𝑖𝑡 • ì 𝑑 is the intersection point.
An essential step in ray tracing is finding ray-primitive intersections. Interestingly, this process bears the similarity of spatial join problems. For example, a complex polygon overlay query can be decomposed into two sub-problems, LSI and PIP [37], which can then be formulated as two RT-friendly problems.
Casting Rays for LSI. Given a base map 𝑅 and a query map 𝑆, LSI finds all the intersections from 𝑅 ⊲⊳ 𝑖𝑛𝑡𝑒𝑟𝑠𝑒𝑐𝑡 𝑆. This can be achieved by ray tracing. We first construct a BVH from an array of AABBs, where each AABB encloses a line segment from 𝑅.
Hit Point AABB line segment from 𝑆 can be represented by a ray. Subsequently, we cast rays generated from 𝑆 and collect all intersections (Figure 4 (b)). The following ray parameters 𝑂, ì 𝑑, 𝑡 𝑚𝑖𝑛 , and 𝑡 𝑚𝑎𝑥 defines a line segment 𝑠 ∈ 𝑆, where 𝑝 1 and 𝑝 2 are the two endpoints of 𝑠.
, which is another endpoint. Thus, we can see that any intersection point 𝑝 𝑥 on 𝑠 can be represented by
, where 𝑡 𝑚𝑖𝑛 ≤ 𝑡 ≤ 𝑡 𝑚𝑎𝑥 , and the ray-AABB intersection test identifies these hit points. Algorithm 1 illustrates the above idea, which is orchestrated by OptiX's callbacks. RayGen shader is the entry of the algorithm that computes 𝑂 and ì 𝑑 using the endpoints. At Line 8, 𝑜𝑝𝑡𝑖𝑥𝑇 𝑟𝑎𝑐𝑒 shoots a ray carrying the line segment ID (per-ray data) from 𝑆 to find intersections. If the ray hits an AABB, a user-defined function IsIntersection is triggered to verify the intersection and calculate the coordinate of intersection point 𝑝 𝑥 and the ray parameter 𝑡 ℎ𝑖𝑡 . This step is necessary because a ray hitting an AABB does not guarantee the ray intersects the line segment enclosed by the AABB. If 𝑝 𝑥 exists, we call 𝑜𝑝𝑡𝑖𝑥𝑅𝑒𝑝𝑜𝑟𝑡𝐼𝑛𝑡𝑒𝑟𝑠𝑒𝑐𝑡𝑖𝑜𝑛, which reports the intersections and extra information with attributes to OptiX.
Upon the reported intersections, if 𝑡 ℎ𝑖𝑡 falls within the interval [𝑡 𝑚𝑖𝑛 , 𝑡 𝑚𝑎𝑥 ], the AnyHit shader is automatically triggered. Remember that each line segment corresponds to an AABB, and we have verified the ray indeed intersects the line segment enclosed by the AABB. We can now safely put the intersection into result set 𝑄. At line 27, we invoke the 𝑜𝑝𝑡𝑖𝑥𝐼𝑔𝑛𝑜𝑟𝑒𝐼𝑛𝑡𝑒𝑟𝑠𝑒𝑐𝑡𝑖𝑜𝑛 function to ask OptiX to continue searching.
Casting Rays for PIP. To begin with, we want to clarify that the classical PIP test determines whether a point lies inside or outside a polygon. This can be addressed with the crossing number algorithm [60]. However, we are solving a more complex PIP problem here. Given a query point and a group of polygons, we want to know which polygon the point falls in. Specifically, this problem is also called "Point Location" in some literature [56]. This operation is the building block of complex spatial queries like polygon overlay [37]. Figure 5 presents an example polygon alongside its in-memory representation, which will be used to demonstrate how the PIP
Algorithm 1: Formulate LSI as a ray-tracing problem Input : 𝑅 -a set of line segments from base map Input :𝑆 -a set of line segments from query map Input :ℎ -the BVH traversal handle built from 𝑅 Output :𝑄 -results of intersections 1 procedure RayGen // Entrypoint, invoked when the program starts 2 𝑡 𝑚𝑖𝑛 , 𝑡 𝑚𝑎𝑥 = {0, 1} 3 for each 𝑠 in 𝑆 do // For each line segment from 𝑆 4 𝑝 1 , 𝑝 2 = 𝐺𝑒𝑡𝐸𝑛𝑑𝑝𝑜𝑖𝑛𝑡𝑠 (𝑠 ) // Get two endpoints from line segment 𝑠 5 𝑂 = 𝑝 1 // Ray origin 6 ì 𝑑 = 𝑝 2 -𝑝 1 // Ray direction 7 // Cast a ray carrying line segment ID of 𝑠 8 𝑜𝑝𝑡𝑖𝑥𝑇 𝑟𝑎𝑐𝑒 (ℎ, 𝑂, ì 𝑑, 𝑡 𝑚𝑖𝑛 , 𝑡 𝑚𝑎𝑥 , 𝑝𝑎𝑦𝑙𝑜𝑎𝑑𝑠 (𝑖𝑑 𝑠 ) ) 9 end for 10 end procedure 11 12 procedure IsIntersection // Invoked when ray potentially hits an AABB 13 𝑂 = 𝑜𝑝𝑡𝑖𝑥𝐺𝑒𝑡𝑅𝑎𝑦𝑂𝑟𝑖𝑔𝑖𝑛 ( ) // Ray origin 14 ì 𝑑 = 𝑜𝑝𝑡𝑖𝑥𝐺𝑒𝑡𝑅𝑎𝑦𝐷𝑖𝑟𝑒𝑐𝑡𝑖𝑜𝑛 ( ) // Ray direction 15 𝑖𝑑 𝑟 = 𝑜𝑝𝑡𝑖𝑥𝐺𝑒𝑡𝑃𝑟𝑖𝑚𝑖𝑡𝑖𝑣𝑒𝐼𝑛𝑑𝑒𝑥 ( ) // Line segment id from 𝑅 16 𝑖𝑑 𝑠 = 𝑜𝑝𝑡𝑖𝑥𝐺𝑒𝑡𝑃𝑎𝑦𝑙𝑜𝑎𝑑 0 ( ) // Line segment id from 𝑆 17 𝑝 = 𝑅 [𝑖𝑑 𝑟 ] ⊲⊳ 𝑖𝑛𝑡𝑒𝑟𝑠𝑒𝑐𝑡𝑠 𝑆 [𝑖𝑑 𝑠 ] // Calculate intersection point 18 if 𝑝 ≠ 𝜙 then 19 𝑡 ℎ𝑖𝑡 = 𝑝𝑥 .𝑥 -𝑂 .𝑥 ì 𝑑 .𝑥 // Calculate 𝑡 ℎ𝑖𝑡 with x-coordinate 20 𝑜𝑝𝑡𝑖𝑥𝑅𝑒𝑝𝑜𝑟𝑡𝐼𝑛𝑡𝑒𝑟𝑠𝑒𝑐𝑡𝑖𝑜𝑛 (𝑡 ℎ𝑖𝑡 , 𝑎𝑡𝑡𝑟𝑖𝑏𝑢𝑡𝑒𝑠 (𝑖𝑑 𝑟 , 𝑖𝑑 𝑠 , 𝑝 ) ) 21 end if 22 end procedure 23 24 procedure AnyHit // Invoked when intersections being reported 25 𝑖𝑑 𝑟 , 𝑖𝑑 𝑠 , 𝑝 𝑥 = 𝑜𝑝𝑡𝑖𝑥𝐺𝑒𝑡𝐴𝑡𝑡𝑟𝑖𝑏𝑢𝑡𝑒 0...2 ( ) // From optixReportIntersection 26 𝑄 = 𝑄 ∪ (𝑖𝑑 𝑟 , 𝑖𝑑 𝑠 , 𝑝 𝑥 ) // Collect the intersection 27 𝑜𝑝𝑡𝑖𝑥𝐼𝑔𝑛𝑜𝑟𝑒𝐼𝑛𝑡𝑒𝑟𝑠𝑒𝑐𝑡𝑖𝑜𝑛 ( ) // Ask OptiX to continue searching 28 end procedure Ray Closest-Hit p 4 p 5 f 2 f 1 f 0 f 0 p 1 p 2 p 3 p 6 p 7 p 8 p 9 p 10 p 11 Query Point (a) PIP query with ray-tracing (b) Chain representation x y O Chain Left Face Right Face Points C 1 f 0 f 1 p 3 , p 2 , p 1 , p 7 C 2 f 2 f 1 p 7 , p 3 C 3 f 2 f 0 p 3 , p 4 , p 5 , p 6 , p 7 C 4 f 0 f 2 p 8 , p 9 , p 10 , p 11 , p 8 p p x Figure 5: Formulating PIP as a ray tracing problem and the chain representation of the left polygon. Notably, 𝐶 4 is a hole, so it only involves two adjacent faces, resulting in a "closed" point sequence (the first point also appears as the last point). works with ray tracing. The point sequence (𝑝 1 , 𝑝 2 , . . . , 𝑝 7 , 𝑝 1 ) consists of the polygon's outline. The line segment (𝑝 7 , 𝑝 3 ) splits the polygon into two facesfoot_3 𝑓 1 and 𝑓 2 . The diamond-shaped polygon (𝑝 8 , 𝑝 9 , 𝑝 10 , 𝑝 11 , 𝑝 8 ) on the bottom is a hole with face 𝑓 0 . A chain is a polyline that belongs to the same adjacent faces. For instance, the polyline (𝑝 3 , 𝑝 2 , 𝑝 1 , 𝑝 7 ) should be organized as a chain because the line segment (𝑝 3 , 𝑝 2 ), (𝑝 2 , 𝑝 1 ), and (𝑝 1 , 𝑝 7 ) share the same left face 𝑓 0 , and right face 𝑓 1 . The order of points determines the adjacent faces of a polyline. For example, in chain 𝐶 2 , the left and right faces are 𝑓 2 and 𝑓 1 , respectively. If the point sequence is (𝑝 3 , 𝑝 7 ), the left face is 𝑓 1 and the right face is 𝑓 2 . The chain representation saves memory by storing common boundaries only once and facilitates the PIP algorithm by providing neighborhood information.
The idea of PIP is to shoot a ray from the query point vertically upwards and to identify the closest line segment hit by the ray [14,37]. This method works both on convex and non-convex polygons, even polygons with holes. The simulation of ray casting is traditionally achieved with grids or trees, which have many shortcomings, as mentioned in §2. Since the PIP problem is essentially a ray tracing problem, solving it with RT Cores is natural. The algorithm works as follows. Figure 5 (a) shows a query 𝑝, locating which polygon the query falls in. We shoot a vertical ray from 𝑝 upwards to identify the closest line segment hit by the ray, which is (𝑝 6 , 𝑝 7 ) in chain 𝐶 3 . Then, we use the chain table (Figure 5 (b)) to determine the face where the query point is located. Since the x-coordinate of point 𝑝 6 is greater than that of 𝑝 7 , we take the left face 𝑓 2 of chain 𝐶 3 , as the polygon ID.
Algorithm 2 describes the implementation of PIP with OptiX. Again, the 𝑅𝑎𝑦𝐺𝑒𝑛 shader is the entry point of the algorithm, starting from casting vertical rays from every point 𝑝 ∈ 𝑆. In IsIntersection shader, we calculate the distance 𝑡 ℎ𝑖𝑡 between the query point 𝑝 and the hit point 𝑝 𝑥 , which is calculated by solving the common solution of line equations (Line 15). Since the ray may hit multiple line segments, we report each 𝑡 ℎ𝑖𝑡 to OptiX by invoking optixRe-portIntersection, asking OptiX to find out the closest hit. When a ray hits the closest line segment, the 𝐶𝑙𝑜𝑠𝑒𝑠𝑡𝐻𝑖𝑡 shader will be called, where we determine which polygon the query falls in. optixGetPayload tells the ID of the query point 𝑖𝑑 𝑝 , and optixGetPrimitiveIndex returns the closest line segment ID 𝑖𝑑 𝑙 . Finally, we retrieve the line segment 𝑙 and query point 𝑝 by the IDs. To tell which polygon 𝑝 falls in, we check the ray casting from the left or right face of 𝑙.
Ray-tracing frameworks like OptiX typically use FP32 for optimal performance [28]. In real-time rendering, such as in video games, the difference in visual quality between lower and higher precision is generally imperceptible to the human eye. However, precision levels are critical in areas like spatial databases. For example, FP32 coordinates offer meter-level distance precision, while FP64 coordinates achieve nanometer-level precision [25], crucial for scientific simulations and engineering analysis. Without improving the precision, the applicability of ray-tracing technology may be severely limited. Precision challenges in graphics, such as self-intersection, and robustness of ray-primitive intersection test, have been wellstudied [22,39]. However, the issue we address here stems from floating-point downcasting and differs from these well-known problems. In this section, we present a simple yet effective method to achieve exact query results for FP64 or higher precision data sources within the constraints of limited hardware support.
To tackle the precision challenges, it is crucial to understand the source of imprecision. Figure 6 represents floating-point numbers Algorithm 2: Formulate PIP as a Ray-Tracing Problem Input : 𝑅 -a set of line segments from base map Input :𝑆 -a set of points from the query map Input :ℎ -the BVH traversal handle built from R Output :𝑄 -Polygon IDs; 𝑄 [𝑖 ] contains the polygon ID of query 𝑆 [𝑖 ] 1 procedure RayGen // Entrypoint, invoked when program starts 2 𝑡 𝑚𝑖𝑛 , 𝑡 𝑚𝑎𝑥 = {0, ∞} 3 for each 𝑝 in 𝑆 do 4 𝑂 = 𝑝 // Query point as origin 5 ì 𝑑.𝑥 𝑦𝑧 = {0, 1, 0} // Towards upward of y-axis 6 𝑜𝑝𝑡𝑖𝑥𝑇 𝑟𝑎𝑐𝑒 (ℎ, 𝑂, ì 𝑑, 𝑡 𝑚𝑖𝑛 , 𝑡 𝑚𝑎𝑥 , 𝑝𝑎𝑦𝑙𝑜𝑎𝑑𝑠 (𝑖𝑑 𝑝 ) ) // Shoot a ray, which carries the point id from 𝑆 7 end for 8 end procedure 9 10 procedure IsIntersection // Invoked when ray potentially hits an AABB 11 𝑖𝑑 𝑙 = 𝑜𝑝𝑡𝑖𝑥𝐺𝑒𝑡𝑃𝑟𝑖𝑚𝑖𝑡𝑖𝑣𝑒𝐼𝑛𝑑𝑒𝑥 ( ) // Line segment id from 𝑅 12 𝑖𝑑 𝑝 = 𝑜𝑝𝑡𝑖𝑥𝐺𝑒𝑡𝑃𝑎𝑦𝑙𝑜𝑎𝑑 0 ( ) // Point id from 𝑆 13 𝑝 = 𝑆 [𝑖𝑑 𝑝 ] 14 𝑙 = 𝑚𝑎𝑘𝑒_𝑣𝑒𝑟𝑡𝑖𝑐𝑎𝑙_𝑙𝑖𝑛𝑒 (𝑝 ) // Generate a vertical line segment crossing 𝑝 15 𝑝 𝑥 = 𝑅 [𝑖𝑑 𝑙 ] ⊲⊳ 𝑖𝑛𝑡𝑒𝑟𝑠𝑒𝑐𝑡𝑠 𝑙 // Intersection point of the line segment from R and 𝑙 16 if 𝑝 𝑥 ≠ 𝜙 then 17 𝑡 ℎ𝑖𝑡 = 𝑝 𝑥 .𝑦 -𝑝.𝑦 // Distance between the ray origin and the hit point 18 𝑜𝑝𝑡𝑖𝑥𝑅𝑒𝑝𝑜𝑟𝑡𝐼𝑛𝑡𝑒𝑟𝑠𝑒𝑐𝑡𝑖𝑜𝑛 (𝑡 ℎ𝑖𝑡 ) // Report intersection to OptiX 19 end if 20 end procedure 21 22 procedure ClosestHit // Invoked when a ray hit the closest line segment 23 𝑖𝑑 𝑙 = 𝑜𝑝𝑡𝑖𝑥𝐺𝑒𝑡𝑃𝑟𝑖𝑚𝑖𝑡𝑖𝑣𝑒𝐼𝑛𝑑𝑒𝑥 ( ) // Line segment id from 𝑅 24 𝑖𝑑 𝑝 = 𝑜𝑝𝑡𝑖𝑥𝐺𝑒𝑡𝑃𝑎𝑦𝑙𝑜𝑎𝑑 0 ( ) // Query point id from 𝑆 25 𝑙 = 𝑅 [𝑖𝑑 𝑙 ] 26 𝑝 1 , 𝑝 2 = 𝐺𝑒𝑡𝐸𝑛𝑑𝑝𝑜𝑖𝑛𝑡𝑠 (𝑙 ) 27 if 𝑝 1 .𝑥 < 𝑝 2 .𝑥 then // Check the line segment direction 28 𝑄 [𝑖𝑑 𝑝 ] = 𝐺𝑒𝑡𝑅𝑖𝑔ℎ𝑡 𝐹𝑎𝑐𝑒 (𝑙 ) 29 else 30 𝑄 [𝑖𝑑 𝑝 ] = 𝐺𝑒𝑡𝐿𝑒 𝑓 𝑡 𝐹𝑎𝑐𝑒 (𝑙 ) 31
end if 32 end procedure 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 1 1 1 1 1 1 0 1 1 0 1 1 0 1 0 0 0 0 1 0 0 1 0 0 1 0 0 0 sign exponent mantissa sign exponent mantissa fp16 fp32 left-over bits Figure 6: Downcasting a 32-bit floating-point to 16-bit. in binary format 6 , consisting of three components: sign, exponent, and [1]. When downcasting, the sign and exponent bits are preserved, ensuring no loss of information. However, the mantissa undergoes truncation. Based on the default rounding rules, known as "roundTiesToEven" [1], the mantissa bits of the low-precision type are left as-is if the discarded bits are below 10 . . . 0, while they are incremented by one if the discarded bits are above 10 . . . 0. In the case of the discarded bits being exactly 10 • • • 0, one of the rounding rules is selected that makes the least significant bit of the mantissa bits zero. Consequently, any downcasted floating-point number can have a value that is less than, greater than, or equal to the value represented by the higher precision type. The process of downcasting high-precision data sources to lowprecision primitives and rays introduces geometric inconsistencies due to the loss of mantissa bits. Figure 7 shows cases that a line segment from dataset 𝑅 intersects two rays casting from dataset 𝑆, resulting in two intersections at x-coordinate 𝑥 1 and 𝑥 2 where 0 < 𝑥 1 < 𝑥 2 . For the sake of simplicity, let's assume that the precision issue only affects the x-axis when constructing the AABB. In Case 1, 𝑥 1 is rounded up, causing the left boundary of the AABB to fail to enclose the left endpoint of the line segment, leading to a false negative outcome. In Case 2, rounding down of 𝑥 2 causes 𝑟𝑎𝑦 2 to miss the AABB, as the right endpoint falls outside the box. In Case 3, neither 𝑟𝑎𝑦 1 nor 𝑟𝑎𝑦 2 intersects the bounding box, resulting in the failure to detect any intersections. These examples highlight the hidden issues that can arise when employing ray-tracing techniques in GIS applications.
1 0 0 1 0 0 1 0 0 0 0 1 1 1 1 1 1 0 1 1 0 1 1 1 0 0 1 0 0 1 0 0 0 An intuitive solution is to create slightly larger AABBs, which can accommodate the line segment after downcasting. However, determining the appropriate size of the enlarged AABB is not trivial. Simply enlarging the AABBs by a small epsilon may still lead to intersections missing, while excessively enlarging them can result in many false positive hits, hurting performance. Formally, let 𝑥 𝑙 denote the x-coordinate of the line segment's left endpoint in high precision. We need to find the maximum value expressible in FP32 as the downcasted left AABB boundary 𝑥 ′ 𝑙 , satisfying 𝑥 ′ 𝑙 ≤ 𝑥 𝑙 . Similarly, 𝑥 𝑟 represents the right endpoint in high precision. We aim to find the minimum value 𝑥 ′ 𝑟 as the downcasted right boundary, ensuring 𝑥 ′ 𝑟 ≥ 𝑥 𝑟 . The same consideration applies to the y-axis. The solution is to create the low-precision AABB that is just barely large enough to fully enclose the high-precision line segment by tweaking mantissa bits. In the scenario depicted in Figure 8, 𝑥 𝑟 is the x-coordinate of the line segment's right endpoint in high precision. Upon downcasting, this value is rounded down, resulting in 𝑥 𝑑 𝑟 , which is less than 𝑥 𝑟 . To determine the smallest possible value 𝑥 ′ 𝑟 as the new right boundary, we only need to add 1 to the mantissa bits after downcasting. This ensures that 𝑥 ′ 𝑟 just exceeds 𝑥 𝑟 . Additionally, to compensate for errors caused by downcasting of the ray origin, we add 1 again to the mantissa bits. To sum up, we add two to the mantissa bits of the right boundary and subtract two from the mantissa bits of the left boundary. Adjusting the mantissa bits may not always be necessary, such as if the 𝑥 𝑟 is rounding up. However, we perform this operation unconditionally to avoid checking all the rounding cases, minimize branches, and enhance execution efficiency on GPUs. Hence, we refer to this approach as Conservative Representation (CR).
The number of AABBs used to build the BVH equals the line segments in the dataset because each line segment is bounded by an individual AABB ( §3.1). This bounding strategy introduces a long construction time and huge memory consumption 7 because the buildup cost is linearly correlated to the number of primitives [83]. One may impulsively use the Ramer-Douglas-Peucker algorithm to simplify polylines and then create an AABB for each simplified polyline [54]. Yet, this algorithm's recursive and sequential nature renders it inefficient for GPU execution. Therefore, we introduce a GPU-optimized method, termed Adaptive Grouping (AG), to group spatially close line segments into the same AABB.
We first introduce the grouping strategy. For two AABBs 𝐵 𝑖 and 𝐵 𝑗 , the union operator ∪ is defined as 𝐵 = 𝐵 𝑖 ∪ 𝐵 𝑗 , where 𝐵 is their MBR, and |𝐵| represents its area. Considering 𝑛 line segments 𝐿 1 , 𝐿 2 , . . . , 𝐿 𝑛 each enclosed by AABBs 𝐵 1 , 𝐵 2 , . . . , 𝐵 𝑛 , we only merge adjacent AABBs in the sequence, such as 𝐵 𝑖 , 𝐵 𝑖+1 , . . . , 𝐵 𝑗 (where 1 ≤ 𝑖 < 𝑗 ≤ 𝑛). This facilitates representing a line segment subset with a simple range 𝑅 [𝑖, 𝑗] instead of explicitly storing each segment's ID. With the range 𝑅, we linearly scan the line segments bounded by 𝐵 to find which line segment is hit by the ray. Merging 𝐵 𝑖 and 𝐵 𝑗 occurs only if the area expansion ratio
where 𝑠 is the merging threshold. We will discuss how the parameter 𝑠 impacts performance and how to find an optimal 𝑠 in §5.5. This approach is inspired by R-Tree's heuristic insertion algorithm [17], which effectively curtails the occurrence of false positives by limiting AABB's dead space.
Figure 9 demonstrates AG's execution. Initially, each line segment (𝐿 0 to 𝐿 7 ) is in a unique group (𝐺 0 to 𝐺 7 ), with a merging threshold 𝑠 = 1.5. In the first iteration (Figure 9 (b)), 𝐿 0 and 𝐿 1 form 𝐺 0 as their area expansion ratio (1.4) is below 𝑠. However, 𝐿 2 and 𝐿 3 remain separate since
The second iteration includes 𝐿 2 in 𝐺 0 and 𝐿 6 in 𝐺 2 , but excludes 𝐿 7 from 𝐺 2 due to a ratio of 1.8, above the threshold. The final iteration merges 𝐿 3 into 𝐺 0 , resulting in three groups: 𝐺 0 with 𝐵 0 and range 𝑅 0 [0, 3], 𝐺 1 with 𝐵 1 and range 𝑅 1 [4,6], and 𝐺 2 with 𝐵 2 and range 𝑅 2 [7,7].
AG can be efficiently implemented on the GPU. We divide 𝐿 segments into ⌈ 𝐿 𝑏 ⌉ partitions, where 𝑏 is the thread block size. 7 One may argue that building the BVH is a one-off cost. A fast buildup allows users to work on exploratory and interactively workloads, e.g., try out different datasets and run some ad-hoc queries quickly.
TB 1 TB 2 (a) Initial Line Segments (b) Iteration 1 (c) Iteration 2
2 TB 1 TB 2 TB 1 TB 2 Figure 9: Execution process of Adaptive Grouping. Initially, a thread block (TB) processes a fixed number of line segments. Spatially close line segments are iteratively put together. Each thread block processes partition 𝑖 containing line segments 𝐿 𝑖 * 𝑏 , 𝐿 𝑖 * 𝑏+1 , . . . , 𝐿 (𝑖+1) * 𝑏 . In the first iteration, thread 𝑗 in thread block 𝑖 attempts to merge two adjacent AABBs 𝐵 𝑖 * 𝑏+𝑗 and 𝐵 𝑖 * 𝑏+𝑗+1 , which encloses line segments 𝐿 𝑖 * 𝑏+𝑗 and 𝐿 𝑖 * 𝑏+𝑗+1 respectively. The thread either produces a merged AABB 𝐵 𝑖 * 𝑏+𝑗 ∪ 𝐵 𝑖 * 𝑏+𝑗+1 and a range 𝑅 [𝑖 * 𝑏 + 𝑗, 𝑖 * 𝑏 + 𝑗 + 1] if the area expanding ratio is within the threshold 𝑠. Otherwise, the AABBs remain unaltered. Subsequent iterations merge two adjacent AABBs from the previous iteration, continuing until no further merges are possible.
Data Storage. The dataset is stored in a planar graph format on host memory and is subsequently decomposed into line segments on device memory for further processing. Each line segment references endpoint IDs instead of storing the coordinates to conserve memory and improve the locality. We pre-compute line equation coefficients from the endpoints. This pre-calculation enables the LSI and PIP queries to access coefficients directly, avoiding redundant calculations. The line segments are then fed into the AG to create AABBs and ranges, followed by applying CR to determine the new boundaries of the AABBs, ensuring precise results.
Coordinates Representation. We employ the fixed-point notation to resolve the accumulated errors from the floating-point arithmetic operations. Given that GPUs do not possess built-in support for fixed-point arithmetic, we convert the floating-point numbers to integers through scaling. This allows all subsequent operations to be performed exclusively on integers, effectively circumventing the issue of precision loss. Scaling to integers also benefits performance because NVIDIA RTX GPUs tend to have very limited FP64 units. Furthermore, we store the intersection points with rational numbers produced by LSI, thereby eliminating the need for fractional coordinates.
Degenerate Cases Handling. In large-scale datasets, degenerate cases such as overlapped line segments for LSI or points on line segments for PIP may occur. These cases must be resolved for correctness. A well-known method for handling degenerate cases is Simulation of Simplicity (SOS) [9]. SOS eliminates all degeneracies by simulating the addition of an infinitesimal value to the coordinates. We employ SOS to handle degenerate cases.
Eliminating Invocation Overhead. Algorithm 1 and 2 conceptually introduce the implementation of RayJoin, where we utilize AnyHit and ClosestHit shaders to collect results. However, employing these shaders may introduce invocation overhead. To mitigate this, both AnyHit and ClosestHit shaders are disabled. We combine the logic of intersection testing and result collection into a single 𝐼𝑠𝐼𝑛𝑡𝑒𝑟𝑠𝑒𝑐𝑡𝑖𝑜𝑛 shader, thereby eliminating the overhead [46].
We also discovered a subset of 𝑆𝑇 functions in geodatabases that may benefit from RT Cores [21]. We only present the ideas of translating them into ray tracing problems. Implementing these operators with RayJoin in geo-database systems is beyond the scope of this paper, and we leave it as future work.
ST_Touches. Given two set of geometric objects 𝐴 and 𝐵, the function 𝑆𝑇 _𝑇𝑜𝑢𝑐ℎ𝑒𝑠 returns true if 𝐶𝑜𝑛𝑑𝑖𝑡𝑖𝑜𝑛 1: 𝐴 and 𝐵 have at least on common point, and 𝐶𝑜𝑛𝑑𝑖𝑡𝑖𝑜𝑛 2: the common points cannot lie in the interior of 𝐴 and 𝐵. Supposing 𝐴 are polygons, we build a BVH from 𝐴. If 𝐵 are points, we cast very short rays as approximations of the points to search hits. If 𝐵 are polygons or polylines, we use Algorithm 1. This fetches all the potential touches, which can be filtered respected with Condition 1. The Condition 2 can be verified with Algorithm 2.
ST_Intersects. This function returns true if the intersection of geometry 𝐴 and 𝐵 is not an empty set. We illustrate the idea due to limited space: a complex geometric object can always be split into primitives (point, line segment). We build a BVH with the primitives decomposed from 𝐴. Then we shoot rays depending on the type of primitives from 𝐵 (a very short ray as a point or a long ray as a line segment). This allows us to collect potentially intersected geometries. Then, we filter the results that satisfy the intersection condition.
ST_Intersection. This function returns the intersected geometry from set 𝐴 and 𝐵. This is one of the most complicated spatial functions and is very compute-intensive. Yet, it benefits from RayJoin. Assuming 𝐴 and 𝐵 are both polygons, which is the most complicated case, we (1) build two BVHs from 𝐴 and 𝐵, (2) run the Algorithm 1 to find all intersections, (3) run Algorithm 2 to determine which polygon the points fall in, and (4) create the resulting polygon by connecting the points of intersection and the points that lie in the interior of polygons. This process is exactly the polygon overlay analysis, which is evaluated in §5.7.
Baselines. Table 1 lists the artifacts we have evaluated. We implemented the uniform grid and LBVH-based solutions on the GPU, which support all the queries [26]. cuSpatial is a GPU-accelerated geospatial analytical platform from NVIDIA [2]. RasterJoin is a raster-based solution on the GPU for spatial aggregation [79]. PSSL is a state-of-the-art plane-sweep-based algorithm for LSI on GPUs [15]. To the best of our knowledge, GLIN is the only learned spatial index that supports indexing geometries with extends [68]. With OpenMP, we use GLIN to implement a parallel LSI query. EPUG-Overlay is a parallel polygon overlay program with a two-level uniform grid [37]. Kinetica and PostGIS are geospatial databases [20]. RasterIntervals uses a raster-based method to accelerate polygon overlay query [16].
Datasets. We collected a diverse set of real-world datasets shown in Table 2, including US-based maps such as County, Zipcode, Block, and Water from ArcGIS Hub [11] and global maps, such as Lakes and Parks from OpenStreetMap (OSM) [10]. We partitioned the OSM datasets based on continents/regions. The spatial join operations were exclusively performed within datasets originating from the same geographic region. 𝑅 ⊲⊳ 𝑆 represents the join operation, where 𝑅 a base map and 𝑆 is a query map. The spatial index is built from 𝑅, and queries are generated from 𝑆. For the scalability experiments, we use Spider to generate synthetic datasets [27]. The polygon size and segment length are set to 0.001 and 10, respectively. The Gaussian distribution is set to 𝜇 = 0.5, 𝜎 = 0.1.
Parameters. Grid-based solutions need a predefined resolution, and trying all the possible resolutions to get the best is prohibitively expensive. By trials, we found that 15000x15000 for the uniform grid yields low running time for all datasets. For EPUG-Overlay, the first-level grid is set to 4000x4000, and the second-level grid is set to 40x40, a typical value used in [37]. RayJoin only needs one parameter, the merging threshold 𝑠, in the adaptive grouping algorithm, which is fixed to 3.5. We discuss how to find an optimal parameter in §5.5.
Environments. RayJoin is implemented in OptiX 8.0 and CUDA 12.3. We evaluate GPU-based solutions on a workstation with Intel Core i9-12900, NVIDIA RTX 3090, and 64 GB memory. cuSpatial version is 23.12.00. For the CPU-based solutions, we employed them on AWS EC2 using an r5.4xlarge instance with 16 vCPUs and 128 GB of RAM. The version of PostGIS is PostGIS-3.3, and Kinetica is v7.1.9.10-ga1. Note that Kinetica has GPU acceleration but it does not support our GPU model and requires a license purchase for self-hosted installation, so we have to evaluate its CPU version on AWS Marketplace Platform.
Table 3 top shows the running time of LSI queries. For the conventional spatial indexes, LBVH is consistently better than the uniform grid. Our analysis attributes the grid-based method's inferior performance to skew distribution issues. For example, in the 𝐶𝑜𝑢𝑛𝑡𝑦 ⊲⊳ 𝑍𝑖𝑝𝑐𝑜𝑑𝑒, almost 100% of cells have workloads (number of intersection tests) under 223K per cell, accounting for 70.97% of total workloads. Yet, about 30% of the workloads are concentrated in just 90 cells (the total number of cells is 225M). PSSL's performance sits between the grid-based and tree-based methods, with its limited parallelism likely hindering the best utilization of the GPU. Notably, GLIN emerges as the least efficient method, with the majority of its processing time devoted to the query stage.
RayJoin demonstrates superior performance across all baselines. Focusing solely on processing time, RayJoin achieves speedups 10.
The bottom half of Table 3 reveals that the uniform grid surpasses LBVH for large datasets, although its performance is hindered by skewed data. The grid's advantage lies in its ability to partition space, allowing for the termination of searches as soon as the closest line segment to a query point is located. In contrast, LBVH requires examining all intersecting line segments to identify the closest one.
RasterJoin, however, shows underperformance compared to both grid and LBVH and runs time out on 𝐿𝐾𝑁 𝐴 ⊲⊳ 𝑃𝐾𝑁 𝐴 in the polygon triangulation stage. For instance, in the 𝐶𝑜𝑢𝑛𝑡𝑦 ⊲⊳ 𝑍𝑖𝑝𝑐𝑜𝑑𝑒 dataset, RasterJoin spent approximately 6878ms on point rendering and 5371ms on polygon rendering, a likely consequence of rendering excessive primitives in large datasets. Notably, NVIDIA's spatial library fell short of the LBVH and ran out-of-memory (OOM) on the 𝐿𝐾𝐸𝑈 ⊲⊳ 𝑃𝐾𝐸𝑈 dataset. Co un ty Zi pc od e Bl oc k W at er LK AF PK AF LK AS PK AS LK AU PK AU LK EU PK EU LK NA PK NA LK SA PK SA 10 1 10 0 10 1 10 2 10 3 10 4 10 5 Speedup (log) OOM Timeout Uniform Grid LBVH RasterJoin cuSpatial PIP is an RT-favorable workload because only the closest hit needs to be found. This property allows us to leverage the ClosestHit shader supported by the ray-tracing hardware to reduce the amount of work instead of traversing all the geometries. If considering the processing time only (ray-tracing), RayJoin exhibited from 4.8x to 74.2x speedups over the best of baselines. As shown in Figure 11, by considering the total running time, RayJoin achieved from 2.1x to 22.2x speeds over the best performance of baselines.
The conservative representation ensures results identical to those of the program using FP64 or higher. To validate its precision, we compared the results of RayJoin to the grid implementation in FP64 on LKNA ⊲⊳ PKNA. The LSI error rate is defined by
, and the PIP error rate is defined by
. Executing LSI and PIP queries with OptiX's default FP32 precision resulted in 745 intersections missing (out of 1,251,343) for LSI and 524 incorrect responses (out of 69,273,661) for PIP, yielding error rates of 0.0595% and 0.0008%, respectively. When using the conservative representation to create the AABBs, both LSI and PIP produce
completely the same answers compared to the FP64 program, and we have verified the correctness of all datasets. These experiments confirm the analysis in §3.2, and our solution eliminates all the false negative cases due to the precision loss.
Parameter Tuning Figure 12 (a) illustrates the variation in BVH construction time, query time, and AG overhead when joining the LKNA with PKNA dataset (Table 2) as parameter 𝑠 changes. When 𝑠 is set to 1, AG is effectively disabled because two line segments are allowed to be put together only if the area of the merged bounding box does not grow, leading to high query performance but also high BVH construction time and memory usage -166ms and 2403MB, respectively (Figure 12 (b-c)). Here, LSI and PIP queries take just 14ms and 18ms.
As 𝑠 increases, construction costs drop significantly, but query times for both LSI and PIP begin to rise. At 𝑠 = 3.5, AG achieves a 29% reduction in total execution time (147ms) and drastically lowers memory usage to about 408MB, which is just 17% of the usage without AG. Beyond 𝑠 = 3.5, total time starts to increase, reaching 155ms at 𝑠 = 5.0. This rise is attributed to AG merging spatially distant bounding boxes, causing a higher false positive rate and more intersection tests.
Figure 12 (b) demonstrates a linear correlation between construction and query times with 𝑠, particularly in the range of 1.5 ≤ 𝑠 ≤ 4. This trend is consistent across various datasets. Optimal parameter selection could be determined through linear regression modeling of the merging threshold-execution time relationship, though this process is omitted here due to space constraints.
We generated synthetic datasets with uniform and Gaussian distributions. The BVH is built from dataset 𝑅 with 5M polygons. A series of datasets 𝑆, ranging from 1M to 5M polygons, were generated to evaluate RayJoin's scalability. To better understand the scalability of RayJoin, adaptive grouping was disabled, and BVH construction time was not included; only query time was reported.
1M 2M 3M 4M 5M Number of Polygons 0 5 10 15 Query Time (ms) (a) Uniform -Query Time LSI -Query Time 1M 2M 3M 4M 5M Number of Polygons 0 10 20 30 40 Query Time (ms) (b) Gaussian -Query Time LSI -Query Time 1M 2M 3M 4M 5M Number of Polygons 0 10 20 30 40 M Intersects/s (c) Uniform -Throughput LSI -Throughput 1M 2M 3M 4M 5M Number of Polygons 0 25 50 75 100 125 M Intersects/s (d) Gaussian -Throughput LSI -Throughput Figure 13 (a)-(d) shows LSI query time and throughput for increasing polygon numbers. Under uniform distribution, 1M polygons result in a 2.9ms query time and 62K intersections. At 3M polygons, the time rises to 6.9ms with 185K intersections, indicating RayJoin's linear scalability, aligning well with the time complexity of 𝑂 (|𝑆 |𝑙𝑜𝑔|𝑅|). Gaussian distributions yield similar patterns but with longer times due to more intersections. When the number of queries increases from 1M to 2M polygons, LSI throughput stabilizes at 29M intersections/s. In the Gaussian dataset, throughput peaks at 99M intersections/s.
For 1M uniform polygons, 8M PIP queries are issued, taking about 5.9ms (Figure 14 (a)). At 2M polygons and 16M queries, the time increases to 12.1ms. In Gaussian distributions, query times are slightly longer due to more ray-primitive intersections. RayJoin
1M 2M 3M 4M 5M Number of Polygons 0 10 20 30 40 Query Time (ms) (a) Uniform -Query Time PIP -Query Time 1M 2M 3M 4M 5M Number of Polygons 0 10 20 30 40 50 Query Time (ms) (b) Gaussian -Query Time PIP -Query Time 1M 2M 3M 4M 5M Number of Polygons 0.0 0.5 1.0 1.5 G Points/s (c) Uniform -Throughput PIP -Throughput 1M 2M 3M 4M 5M Number of Polygons 0.0 0.5 1.0 1.5 2.0 G Points/s (d) Gaussian -Throughput PIP -Throughput shows linear scalability with query numbers. Figure 14 (c), (d) indicates a peak PIP throughput of around 1.5G points/s for both distributions, affirming RayJoin is highly scalable as throughput remains consistent despite increased queries.
The polygon overlay query identifies and outputs the overlapped area of intersected polygons from two maps, which can be achieved by combining the LSI and PIP query results ( §2.1). The overlay processing is extremely time-consuming due to the complexity of the overlay algorithm and the irregularity of the execution.
Co un ty Zi pc od e Bl oc k W at er LK AF PK AF LK AS PK AS LK AU PK AU LK EU PK EU LK NA PK NA LK SA PK SA 10 1 10 0 10 1 10 2 10 3 10 4 10 5 Speedup (log) OOM PostGIS Kinetica EPUG-Overlay RasterIntervals Uniform Grid LBVH Table 4 shows the polygon overlay running time. It is evident that none of the existing solutions other than RayJoin take under one second for all datasets, even the highly optimized GPU solutions (uniform grid and LBVH). The longest execution time of PostGIS is about 234s on the 𝐵𝑙𝑜𝑐𝑘 ⊲⊳ 𝑊 𝑎𝑡𝑒𝑟 dataset. Kinetica is faster than PostGIS because it is memory-based, only taking about 53s. EPUG-Overlay typically takes more than hundreds of seconds for fairly large datasets. Although RasterIntervals is fast at the processing stage, it can spend thousands of seconds on preprocessing. The uniform grid takes up to 9.4s for processing 𝐿𝐾𝐸𝑈 ⊲⊳ 𝑃𝐾𝐸𝑈 , and LBVH even runs OOM on this dataset. LBVH takes up to 58.4s to process the 𝐿𝐾𝑁 𝐴 ⊲⊳ 𝑃𝐾𝑁 𝐴 dataset due to irregular random memory access patterns of tree-based structure and the inability to terminate the PIP search early.
RayJoin outperformed all the state-of-the-art solutions on both CPUs and GPUs. By leveraging the RT Cores, RayJoin efficiently transforms the polygon overlay problem into two RT-friendly problems. The RT Cores overcome the long memory access latency and irregular memory access patterns while benefiting from the logarithmic complexity of tree-based indexes, resulting in extremely low query time. As summarized in Figure 15, RayJoin achieved from 3.0x to 28.3x speed up over the best of the baselines.
Spatial Indexing on Advanced Hardware. Geometric Performance Primitives (GPP) is a robust polygon overlay library on GPUs with a uniform grid [6]. Liu et al. proposed a filtering technique to accelerate spatial join with fine-grained tiles [33,34]. Rayhan et al. designed an R-tree that leverages SIMD instructions [55]. SwiftSpatial is a spatial join accelerator on FPGAs [23]. Note that we target the same problem but with different hardware. In short, RayJoin is an index-based solution uniquely exploiting ray-tracing hardware's capability to deliver real-time spatial join.
Rasterization-based Techniques. RasterJoin is a rasterizationbased method for spatial aggregation by leveraging the rendering pipeline [79]. RasterIntervals is the state-of-the-art polygon overlay method with the approximation technique to reduce the refinement cost [16]. PixelBox employs a rasterization-based method to compute the Jaccard Similarity [70]. IDEAL is a vector-raster hybrid method for PIP and point-to-polygon distance calculation [65]. Compared to rasterization-based methods, RayJoin does not need any preprocessing steps to build up rasterized representations.
Parallel Plane Sweep Algorithms. McKenney et al. proposed a parallel plane sweep algorithm on CPUs [41] that dynamically segments the input to spatial operations and executes them in parallel. PSSL is the plane sweep-based solution for the LSI query on the GPU [15]. The plane sweep algorithms are challenging to parallelize, as demonstrated by our evaluations, which clearly indicate that they exhibit significantly lower performance than other solutions.
Accelerating Non-graphics Workload with RT Cores. An early work on repurposing RT Cores, proposed by Salmon et al., involves accelerating Monte Carlo particle transport simulations [57]. RTNN adopted spatial reordering and query partitioning optimizations for K nearest neighbor (KNN). TrueKNN supports arbitrary radius KNN-search on RT Cores [46]. JUNO is the first work supporting high-dimensional KNN search on RT Cores [35]. A method presented in [38] supports non-Euclidean distance measures for KNN. A method proposed in [30] accelerates the PIP test by transforming polygons into 3D representations, which is different from our method and has accuracy limitations. Accelerating Tet-Mesh Point Location with RT Cores is discussed in [67]. Literature [24] explores using RT Cores to serve as a B-Tree. RTScan leverages RT Cores to accelerate index scans [36]. Accelerating Range Minimum Queries with RT Cores is discussed in [42]. None of these works address the spatial join problem that our paper targets, which presents unique challenges such as problem formulation, precision preservation, and performance optimization.
To understand dataset R, considering the provinces/states (the three simple polygons) share the common boundaries forms into a country (the whole polygon).
AABB is a synonym of MBR in this paper.
We use the word "primitive" interchangeably with "geometry." In the context of BVH, we prefer the word "primitive" because they are conventionally used together in computer graphics.
The term "face" defines a planar region enclosed by the boundaries. The hole and region that is not enclosed are defined as the exterior face 𝑓 0 .
For ease of representation, we demonstrate the process of downcasting from 32 bits to 16 bits instead of 64 bits to 32 bits.