Repository: rapidsai/cuspatial (language: python)

Issue:
[BUG]: Strange behavior with quadtree point in polygon
### Version

23.04

### On which installation method(s) does this occur?

Rapids-Compose

### Describe the issue

quadtree pip seems to experience a silent OOM or other error in certain larger data cases. This particular example requires the availability of our demo datasets `taxi2015.csv` and `taxi_zones.zip`. I'm not sure of another way to reproduce it, so I'm including the full example here. It doesn't always appear so I'm just remembering the issue here.

As you can see from the Relevant log output, quadtree returns 3.7m samples when only the first 120 polygons in `zones` are used. If the entire `zones` polygons list (263 polygons) is used, _something_ happens and only 16314 rows are returned. This is unexpected behavior and needs investigation.

### Minimum reproducible example

```shell
# https://github.com/rapidsai/cuspatial/issues/890
import cuspatial
import cudf
import geopandas
host_zones = geopandas.read_file('taxi_zones.zip')
host_lonlat = host_zones.to_crs(epsg=4326)
zones = cuspatial.from_geopandas(host_lonlat)
taxi2015 = cudf.read_csv('taxi2015.csv')
def quadtree(polygons, points):
    poly_points_x = polygons.polygons.x
    poly_points_y = polygons.polygons.y
    poly_offsets = polygons.polygons.part_offset
    poly_ring_offsets = polygons.polygons.ring_offset
    test_points_x = points.points.x
    test_points_y = points.points.y
    scale = 50
    max_depth = 7
    min_size = 125
    x_max = poly_points_x.max()
    x_min = poly_points_x.min()
    y_max = poly_points_y.max()
    y_min = poly_points_y.min()
    point_indices, quadtree = cuspatial.quadtree_on_points(
        test_points_x,
        test_points_y,
        x_min,
        x_max,
        y_min,
        y_max,
        scale,
        max_depth,
        min_size,
    )
    poly_bboxes = cuspatial.polygon_bounding_boxes(
        poly_offsets, poly_ring_offsets, poly_points_x, poly_points_y
    )
    intersections = cuspatial.join_quadtree_and_bounding_boxes(
        quadtree, poly_bboxes, x_min, x_max, y_min, y_max, scale, max_depth
    )
    polygons_and_points = cuspatial.quadtree_point_in_polygon(
        intersections,
        quadtree,
        point_indices,
        test_points_x,
        test_points_y,
        poly_offsets,
        poly_ring_offsets,
        poly_points_x,
        poly_points_y,
    )
    return polygons_and_points
def make_geoseries_from_lon_lat(lon, lat):
    # Scatter the two columns into one column
    assert len(lon) == len(lat)
    xy = cudf.Series(cp.zeros(len(lon) * 2))
    xy[::2] = lon
    xy[1::2] = lat

    return cuspatial.GeoSeries(cuspatial.core._column.geocolumn.GeoColumn._from_points_xy(xy._column))
dropoffs = make_geoseries_from_lon_lat(
    taxi2015['dropoff_longitude'],
    taxi2015['dropoff_latitude']
)
print(quadtree(zones['geometry'].iloc[0:120], dropoffs))
print(quadtree(zones['geometry'], dropoffs))
```


### Relevant log output

```shell
polygon_index  point_index
0                    0          116
1                    0          387
2                    0          685
3                    0         2607
4                    0         3141
...                ...          ...
3686177            167     12304018
3686178            167     12323531
3686179            167     12351800
3686180            167     12444884
3686181            167     12484251

[3686182 rows x 2 columns]
       polygon_index  point_index
0                  0          116
1                  0          387
2                  0          685
3                  0         2608
4                  0         3142
...              ...          ...
16309              0     12504645
16310              1      7456107
16311              1      7530752
16312              1      7704910
16313              1     11938181

[16314 rows x 2 columns]
```


### Environment details

```shell
rapids-compose 23.04
```


### Other/Misc.

_No response_


Discussion / hints:
If the point coordiantes are also using projection EPSG:4326 (lon/lat), it does not make sense to use a scale of 50. The width/height of NYC bounding boxes is no more than 1 degree in lon/lat. 
The issue I'm asking about has to do with the hidden error that causes the full dataset to return fewer points than taking a partial slice of it.

`print(quadtree(zones['geometry'], dropoffs))` returns only 16313 points
`print(quadtree(zones['geometry'].iloc[0:120], dropoffs))` returns 3686181 points
`print(quadtree(zones['geometry'].iloc[0:160], dropoffs))` returns an OOM error
`print(quadtree(zones['geometry'].iloc[0:len(zones)-9]` also quietly returns a smaller number of points, or any number down to 1. If I use `len(zones)-10` I get an OOM error.

Fix patch:
diff --git a/cpp/benchmarks/CMakeLists.txt b/cpp/benchmarks/CMakeLists.txt
index 99780a677..10e626f3f 100644
--- a/cpp/benchmarks/CMakeLists.txt
+++ b/cpp/benchmarks/CMakeLists.txt
@@ -1,5 +1,5 @@
 #=============================================================================
-# Copyright (c) 2019-2021, NVIDIA CORPORATION.
+# Copyright (c) 2019-2024, NVIDIA CORPORATION.
 #
 # Licensed under the Apache License, Version 2.0 (the "License");
 # you may not use this file except in compliance with the License.
@@ -23,6 +23,17 @@ add_library(cuspatial_benchmark_common OBJECT
 
 target_compile_features(cuspatial_benchmark_common PUBLIC cxx_std_17 cuda_std_17)
 
+set_target_properties(cuspatial_benchmark_common
+    PROPERTIES RUNTIME_OUTPUT_DIRECTORY "$<BUILD_INTERFACE:${CUSPATIAL_BINARY_DIR}/benchmarks>"
+               INSTALL_RPATH "\$ORIGIN/../../../lib"
+               CXX_STANDARD                        17
+               CXX_STANDARD_REQUIRED               ON
+               CUDA_STANDARD                       17
+               CUDA_STANDARD_REQUIRED              ON
+               POSITION_INDEPENDENT_CODE           ON
+               INTERFACE_POSITION_INDEPENDENT_CODE ON
+)
+
 target_link_libraries(cuspatial_benchmark_common
     PUBLIC benchmark::benchmark
            cudf::cudftestutil
@@ -43,6 +54,10 @@ function(ConfigureBench CMAKE_BENCH_NAME)
     set_target_properties(${CMAKE_BENCH_NAME}
         PROPERTIES RUNTIME_OUTPUT_DIRECTORY "$<BUILD_INTERFACE:${CUSPATIAL_BINARY_DIR}/benchmarks>"
                    INSTALL_RPATH "\$ORIGIN/../../../lib"
+                   CXX_STANDARD 17
+                   CXX_STANDARD_REQUIRED ON
+                   CUDA_STANDARD 17
+                   CUDA_STANDARD_REQUIRED ON
         )
     target_link_libraries(${CMAKE_BENCH_NAME} PRIVATE benchmark::benchmark_main cuspatial_benchmark_common)
     install(
@@ -61,7 +76,11 @@ function(ConfigureNVBench CMAKE_BENCH_NAME)
     ${CMAKE_BENCH_NAME}
     PROPERTIES RUNTIME_OUTPUT_DIRECTORY "$<BUILD_INTERFACE:${CUSPATIAL_BINARY_DIR}/benchmarks>"
                INSTALL_RPATH "\$ORIGIN/../../../lib"
-  )
+               CXX_STANDARD 17
+               CXX_STANDARD_REQUIRED ON
+               CUDA_STANDARD 17
+               CUDA_STANDARD_REQUIRED ON
+)
   target_link_libraries(
     ${CMAKE_BENCH_NAME} PRIVATE cuspatial_benchmark_common nvbench::main
   )
diff --git a/cpp/include/cuspatial/detail/join/quadtree_point_in_polygon.cuh b/cpp/include/cuspatial/detail/join/quadtree_point_in_polygon.cuh
index 0b8f23e19..68a2e9596 100644
--- a/cpp/include/cuspatial/detail/join/quadtree_point_in_polygon.cuh
+++ b/cpp/include/cuspatial/detail/join/quadtree_point_in_polygon.cuh
@@ -24,14 +24,16 @@
 #include <cuspatial/range/multipolygon_range.cuh>
 #include <cuspatial/traits.hpp>
 
+#include <rmm/cuda_device.hpp>
 #include <rmm/device_uvector.hpp>
 #include <rmm/exec_policy.hpp>
 #include <rmm/resource_ref.hpp>
 
 #include <thrust/iterator/permutation_iterator.h>
+#include <thrust/iterator/transform_iterator.h>
 #include <thrust/scan.h>
 
-#include <cstdint>
+#include <limits>
 
 namespace cuspatial {
 namespace detail {
@@ -57,7 +59,7 @@ struct compute_poly_and_point_indices {
   using IndexType = iterator_value_type<QuadOffsetsIterator>;
 
   inline thrust::tuple<IndexType, IndexType> __device__
-  operator()(IndexType const global_index) const
+  operator()(std::uint64_t const global_index) const
   {
     auto const [quad_poly_index, local_point_index] =
       get_quad_and_local_point_indices(global_index, point_offsets_begin, point_offsets_end);
@@ -118,16 +120,26 @@ std::pair<rmm::device_uvector<IndexType>, rmm::device_uvector<IndexType>> quadtr
 
   auto num_poly_quad_pairs = std::distance(poly_indices_first, poly_indices_last);
 
-  auto quad_lengths_iter =
-    thrust::make_permutation_iterator(quadtree.length_begin(), quad_indices_first);
+  // The quadtree length is an iterator of uint32_t, but we have to transform into uint64_t values
+  // so the thrust::inclusive_scan accumulates into uint64_t outputs. Changing the output iterator
+  // to uint64_t isn't sufficient to achieve this behavior.
+  auto quad_lengths_iter = thrust::make_transform_iterator(
+    thrust::make_permutation_iterator(quadtree.length_begin(), quad_indices_first),
+    cuda::proclaim_return_type<std::uint64_t>([] __device__(IndexType const& i) -> std::uint64_t {
+      return static_cast<std::uint64_t>(i);
+    }));
 
   auto quad_offsets_iter =
     thrust::make_permutation_iterator(quadtree.offset_begin(), quad_indices_first);
 
-  // Compute a "local" set of zero-based point offsets from number of points in each quadrant
+  // Compute a "local" set of zero-based point offsets from the number of points in each quadrant.
+  //
   // Use `num_poly_quad_pairs + 1` as the length so that the last element produced by
   // `inclusive_scan` is the total number of points to be tested against any polygon.
-  rmm::device_uvector<IndexType> local_point_offsets(num_poly_quad_pairs + 1, stream);
+  //
+  // Accumulate into uint64_t, because the prefix sums can overflow the size of uint32_t
+  // when testing a large number of polygons against a large quadtree.
+  rmm::device_uvector<std::uint64_t> local_point_offsets(num_poly_quad_pairs + 1, stream);
 
   // inclusive scan of quad_lengths_iter
   thrust::inclusive_scan(rmm::exec_policy(stream),
@@ -136,21 +148,27 @@ std::pair<rmm::device_uvector<IndexType>, rmm::device_uvector<IndexType>> quadtr
                          local_point_offsets.begin() + 1);
 
   // Ensure local point offsets starts at 0
-  IndexType init{0};
+  std::uint64_t init{0};
   local_point_offsets.set_element_async(0, init, stream);
 
   // The last element is the total number of points to test against any polygon.
   auto num_total_points = local_point_offsets.back_element(stream);
 
-  // Allocate the output polygon and point index pair vectors
-  rmm::device_uvector<IndexType> poly_indices(num_total_points, stream);
-  rmm::device_uvector<IndexType> point_indices(num_total_points, stream);
-
-  auto poly_and_point_indices =
-    thrust::make_zip_iterator(poly_indices.begin(), point_indices.begin());
-
-  // Enumerate the point X/Ys using the sorted `point_indices` (from quadtree construction)
-  auto point_xys_iter = thrust::make_permutation_iterator(points_first, point_indices_first);
+  // The largest supported input size for thrust::count_if/copy_if is INT32_MAX.
+  // This functor iterates over the input space and processes up to INT32_MAX elements at a time.
+  std::uint64_t max_points_to_test = std::numeric_limits<std::int32_t>::max();
+  auto count_in_chunks             = [&](auto const& func) {
+    std::uint64_t memo{};
+    for (std::uint64_t offset{0}; offset < num_total_points; offset += max_points_to_test) {
+      memo += func(memo, offset, std::min(max_points_to_test, num_total_points - offset));
+    }
+    return memo;
+  };
+
+  detail::test_poly_point_intersection test_poly_point_pair{
+    // Enumerate the point X/Ys using the sorted `point_indices` (from quadtree construction)
+    thrust::make_permutation_iterator(points_first, point_indices_first),
+    polygons};
 
   // Compute the combination of polygon and point index pairs. For each polygon/quadrant pair,
   // enumerate pairs of (poly_index, point_index) for each point in each quadrant.
@@ -163,28 +181,57 @@ std::pair<rmm::device_uvector<IndexType>, rmm::device_uvector<IndexType>> quadtr
   //     pp_pairs.append((polygon, point))
   // ```
   //
-  auto global_to_poly_and_point_indices = detail::make_counting_transform_iterator(
-    0,
-    detail::compute_poly_and_point_indices{quad_offsets_iter,
-                                           local_point_offsets.begin(),
-                                           local_point_offsets.end(),
-                                           poly_indices_first});
-
-  // Compute the number of intersections by removing (poly, point) pairs that don't intersect
-  auto num_intersections = thrust::distance(
-    poly_and_point_indices,
-    thrust::copy_if(rmm::exec_policy(stream),
-                    global_to_poly_and_point_indices,
-                    global_to_poly_and_point_indices + num_total_points,
-                    poly_and_point_indices,
-                    detail::test_poly_point_intersection{point_xys_iter, polygons}));
-
-  poly_indices.resize(num_intersections, stream);
-  poly_indices.shrink_to_fit(stream);
-  point_indices.resize(num_intersections, stream);
-  point_indices.shrink_to_fit(stream);
-
-  return std::pair{std::move(poly_indices), std::move(point_indices)};
+  auto global_to_poly_and_point_indices = [&](auto offset = 0) {
+    return detail::make_counting_transform_iterator(
+      offset,
+      detail::compute_poly_and_point_indices{quad_offsets_iter,
+                                             local_point_offsets.begin(),
+                                             local_point_offsets.end(),
+                                             poly_indices_first});
+  };
+
+  auto run_quadtree_point_in_polygon = [&](auto output_size) {
+    // Allocate the output polygon and point index pair vectors
+    rmm::device_uvector<IndexType> poly_indices(output_size, stream);
+    rmm::device_uvector<IndexType> point_indices(output_size, stream);
+
+    auto num_intersections = count_in_chunks([&](auto memo, auto offset, auto size) {
+      auto poly_and_point_indices =
+        thrust::make_zip_iterator(poly_indices.begin(), point_indices.begin()) + memo;
+      // Remove (poly, point) pairs that don't intersect
+      return thrust::distance(poly_and_point_indices,
+                              thrust::copy_if(rmm::exec_policy(stream),
+                                              global_to_poly_and_point_indices(offset),
+                                              global_to_poly_and_point_indices(offset) + size,
+                                              poly_and_point_indices,
+                                              test_poly_point_pair));
+    });
+
+    if (num_intersections < output_size) {
+      poly_indices.resize(num_intersections, stream);
+      point_indices.resize(num_intersections, stream);
+      poly_indices.shrink_to_fit(stream);
+      point_indices.shrink_to_fit(stream);
+    }
+
+    return std::pair{std::move(poly_indices), std::move(point_indices)};
+  };
+
+  try {
+    // First attempt to run the hit test assuming allocating space for all possible intersections
+    // fits into the available memory.
+    return run_quadtree_point_in_polygon(num_total_points);
+  } catch (rmm::out_of_memory const&) {
+    // If we OOM the first time, pre-compute the number of hits and allocate only that amount of
+    // space for the output buffers. This halves performance, but it should at least return valid
+    // results.
+    return run_quadtree_point_in_polygon(count_in_chunks([&](auto memo, auto offset, auto size) {
+      return thrust::count_if(rmm::exec_policy(stream),
+                              global_to_poly_and_point_indices(offset),
+                              global_to_poly_and_point_indices(offset) + size,
+                              test_poly_point_pair);
+    }));
+  }
 }
 
 }  // namespace cuspatial
diff --git a/cpp/include/cuspatial/detail/range/multilinestring_range.cuh b/cpp/include/cuspatial/detail/range/multilinestring_range.cuh
index b9b53bfc0..03ad0fe27 100644
--- a/cpp/include/cuspatial/detail/range/multilinestring_range.cuh
+++ b/cpp/include/cuspatial/detail/range/multilinestring_range.cuh
@@ -1,5 +1,5 @@
 /*
- * Copyright (c) 2022-2023, NVIDIA CORPORATION.
+ * Copyright (c) 2022-2024, NVIDIA CORPORATION.
  *
  * Licensed under the Apache License, Version 2.0 (the "License");
  * you may not use this file except in compliance with the License.
@@ -74,6 +74,7 @@ template <typename GeometryIterator, typename PartIterator, typename VecIterator
 class multilinestring_range;
 
 template <typename GeometryIterator, typename PartIterator, typename VecIterator>
+CUSPATIAL_HOST_DEVICE
 multilinestring_range<GeometryIterator, PartIterator, VecIterator>::multilinestring_range(
   GeometryIterator geometry_begin,
   GeometryIterator geometry_end,
diff --git a/cpp/include/cuspatial/geometry/box.hpp b/cpp/include/cuspatial/geometry/box.hpp
index 1041c4de2..4a9f97639 100644
--- a/cpp/include/cuspatial/geometry/box.hpp
+++ b/cpp/include/cuspatial/geometry/box.hpp
@@ -1,5 +1,5 @@
 /*
- * Copyright (c) 2022, NVIDIA CORPORATION.
+ * Copyright (c) 2022-2024, NVIDIA CORPORATION.
  *
  * Licensed under the Apache License, Version 2.0 (the "License");
  * you may not use this file except in compliance with the License.
@@ -40,9 +40,9 @@ class alignas(sizeof(Vertex)) box {
 
  private:
   /**
-   * @brief Output stream operator for `vec_2d<T>` for human-readable formatting
+   * @brief Output stream operator for `box<T>` for human-readable formatting
    */
-  friend std::ostream& operator<<(std::ostream& os, cuspatial::box<T> const& b)
+  friend std::ostream& operator<<(std::ostream& os, cuspatial::box<T, Vertex> const& b)
   {
     return os << "{" << b.v1 << ", " << b.v2 << "}";
   }


Test patch:
diff --git a/cpp/include/cuspatial_test/geometry_generator.cuh b/cpp/include/cuspatial_test/geometry_generator.cuh
index a3efa5812..8c17ff291 100644
--- a/cpp/include/cuspatial_test/geometry_generator.cuh
+++ b/cpp/include/cuspatial_test/geometry_generator.cuh
@@ -33,6 +33,8 @@
 
 #include <ranger/ranger.hpp>
 
+#include <cmath>
+
 namespace cuspatial {
 namespace test {
 
@@ -399,8 +401,9 @@ auto generate_multipoint_array(multipoint_generator_parameter<T> params,
                    std::size_t{0},
                    params.num_points_per_multipoints);
 
-  auto engine_x = deterministic_engine(params.num_points());
-  auto engine_y = deterministic_engine(2 * params.num_points());
+  auto golden_ratio = (1 + std::sqrt(T{5})) / 2;
+  auto engine_x     = deterministic_engine(golden_ratio * params.num_points());
+  auto engine_y     = deterministic_engine((1 / golden_ratio) * params.num_points());
 
   auto x_dist = make_uniform_dist(params.lower_left.x, params.upper_right.x);
   auto y_dist = make_uniform_dist(params.lower_left.y, params.upper_right.y);
diff --git a/cpp/include/cuspatial_test/vector_factories.cuh b/cpp/include/cuspatial_test/vector_factories.cuh
index a772cd6cf..9d3e8a31a 100644
--- a/cpp/include/cuspatial_test/vector_factories.cuh
+++ b/cpp/include/cuspatial_test/vector_factories.cuh
@@ -363,6 +363,17 @@ class multipoint_array {
       _geometry_offsets.begin(), _geometry_offsets.end(), _coordinates.begin(), _coordinates.end()};
   }
 
+  /**
+   * @brief Copy the offset arrays to host.
+   */
+  auto to_host() const
+  {
+    auto geometry_offsets   = cuspatial::test::to_host<geometry_t>(_geometry_offsets);
+    auto coordinate_offsets = cuspatial::test::to_host<coord_t>(_coordinates);
+
+    return std::tuple{geometry_offsets, coordinate_offsets};
+  }
+
   /// Release ownership
   auto release() { return std::pair{std::move(_geometry_offsets), std::move(_coordinates)}; }
 
diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt
index 7f7cb112a..f6752c860 100644
--- a/cpp/tests/CMakeLists.txt
+++ b/cpp/tests/CMakeLists.txt
@@ -1,5 +1,5 @@
 ﻿#=============================================================================
-# Copyright (c) 2019-2023, NVIDIA CORPORATION.
+# Copyright (c) 2019-2024, NVIDIA CORPORATION.
 #
 # Licensed under the Apache License, Version 2.0 (the "License");
 # you may not use this file except in compliance with the License.
@@ -158,6 +158,9 @@ ConfigureTest(JOIN_POINT_IN_POLYGON_SMALL_TEST_EXP
 ConfigureTest(JOIN_POINT_IN_POLYGON_LARGE_TEST_EXP
     join/quadtree_point_in_polygon_test_large.cu)
 
+ConfigureTest(JOIN_POINT_IN_POLYGON_OOM_TEST_EXP
+    join/quadtree_point_in_polygon_test_oom.cu)
+
 ConfigureTest(JOIN_POINT_TO_LINESTRING_SMALL_TEST_EXP
     join/quadtree_point_to_nearest_linestring_test_small.cu)
 
diff --git a/cpp/tests/join/quadtree_point_in_polygon_test_oom.cu b/cpp/tests/join/quadtree_point_in_polygon_test_oom.cu
new file mode 100644
index 000000000..2f66918e0
--- /dev/null
+++ b/cpp/tests/join/quadtree_point_in_polygon_test_oom.cu
@@ -0,0 +1,140 @@
+/*
+ * Copyright (c) 2023-2024, NVIDIA CORPORATION.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include <cuspatial_test/base_fixture.hpp>
+#include <cuspatial_test/geometry_generator.cuh>
+#include <cuspatial_test/vector_factories.cuh>
+
+#include <cuspatial/bounding_boxes.cuh>
+#include <cuspatial/detail/algorithm/is_point_in_polygon.cuh>
+#include <cuspatial/geometry/box.hpp>
+#include <cuspatial/geometry/vec_2d.hpp>
+#include <cuspatial/iterator_factory.cuh>
+#include <cuspatial/spatial_join.cuh>
+
+#include <gtest/gtest.h>
+
+/*
+ * The test uses the same quadtree structure as in pip_refine_test_small. However, the number of
+ * randomly generated points under all quadrants (min_size) are increased to be more than the
+ * number of threads per-block.
+ */
+
+template <typename T>
+struct PIPRefineTestLarge : public cuspatial::test::BaseFixture {};
+
+template <typename T>
+struct test_point_in_poly_functor {
+  cuspatial::multipolygon_range<std::size_t*, std::size_t*, std::size_t*, cuspatial::vec_2d<T>*>
+    polys;
+  __device__ inline T operator()(cuspatial::vec_2d<T> point)
+  {
+    auto it = cuspatial::make_counting_transform_iterator(0, [&](auto i) { return polys[i][0]; });
+    return thrust::count_if(thrust::seq, it, it + polys.num_multipolygons(), [&](auto poly) {
+      return cuspatial::is_point_in_polygon(point, poly);
+    });
+  }
+};
+
+using TestTypes = ::testing::Types<float, double>;
+
+TYPED_TEST_CASE(PIPRefineTestLarge, TestTypes);
+
+TYPED_TEST(PIPRefineTestLarge, TestOOM)
+{
+  using T = TypeParam;
+  using cuspatial::vec_2d;
+  using cuspatial::test::make_device_vector;
+
+  using point_t = typename cuspatial::test::multipoint_array<rmm::device_uvector<std::size_t>,
+                                                             rmm::device_uvector<vec_2d<T>>>;
+
+  using polys_t = typename cuspatial::test::multipolygon_array<rmm::device_uvector<std::size_t>,
+                                                               rmm::device_uvector<std::size_t>,
+                                                               rmm::device_uvector<std::size_t>,
+                                                               rmm::device_uvector<vec_2d<T>>>;
+
+  vec_2d<T> v_min{0.0, 0.0};
+  vec_2d<T> v_max{1'000.0, 1'000.0};
+  T const scale{1.0};
+  std::uint8_t const max_depth{15};
+  std::uint32_t const min_size{10'000};
+
+  std::size_t const num_points{1'000'000};
+  std::size_t const num_polys{24'000};
+
+  rmm::device_uvector<vec_2d<T>> points = [&]() {
+    point_t points = cuspatial::test::generate_multipoint_array(
+      cuspatial::test::multipoint_generator_parameter<T>{1, num_points, v_min, v_max},
+      this->stream());
+    return points.release().second;
+  }();
+
+  polys_t multipoly_array = cuspatial::test::generate_multipolygon_array(
+    cuspatial::test::multipolygon_generator_parameter<T>{
+      num_polys,
+      1,
+      0,
+      4,
+      vec_2d<T>{(v_max - v_min).x / 2, (v_max - v_min).y / 2},
+      (v_max - v_min).x / 8},
+    this->stream());
+
+  auto multipolygons = multipoly_array.range();
+
+  auto expected_size = thrust::count_if(rmm::exec_policy(this->stream()),
+                                        points.begin(),
+                                        points.end(),
+                                        test_point_in_poly_functor<T>{multipolygons});
+
+  auto [point_indices, quadtree] = cuspatial::quadtree_on_points(
+    points.begin(), points.end(), v_min, v_max, scale, max_depth, min_size, this->stream());
+
+  auto bboxes =
+    rmm::device_uvector<cuspatial::box<T>>(multipolygons.num_polygons(), this->stream());
+
+  cuspatial::polygon_bounding_boxes(multipolygons.part_offset_begin(),
+                                    multipolygons.part_offset_end(),
+                                    multipolygons.ring_offset_begin(),
+                                    multipolygons.ring_offset_end(),
+                                    multipolygons.point_begin(),
+                                    multipolygons.point_end(),
+                                    bboxes.begin(),
+                                    T{0},
+                                    this->stream());
+
+  EXPECT_GT(bboxes.size(), 0);
+
+  auto [poly_indices, quad_indices] = cuspatial::join_quadtree_and_bounding_boxes(
+    quadtree, bboxes.begin(), bboxes.end(), v_min, scale, max_depth, this->stream());
+
+  auto [actual_poly_indices, actual_point_indices] =
+    cuspatial::quadtree_point_in_polygon(poly_indices.begin(),
+                                         poly_indices.end(),
+                                         quad_indices.begin(),
+                                         quadtree,
+                                         point_indices.begin(),
+                                         point_indices.end(),
+                                         points.begin(),
+                                         multipolygons,
+                                         this->stream());
+
+  EXPECT_GT(actual_point_indices.size(), 0);
+  EXPECT_GT(actual_poly_indices.size(), 0);
+
+  EXPECT_EQ(actual_point_indices.size(), expected_size);
+  EXPECT_EQ(actual_poly_indices.size(), expected_size);
+}