Add find_point_on_segment internal API

cpp/include/cuspatial/detail/utility/floating_point.cuh

@@ -127,25 +127,5 @@ bool CUSPATIAL_HOST_DEVICE float_equal(T const& flhs, T const& frhs)
                                 : (rhsbiased - lhsbiased) <= max_ulp;
 }
 
-/**
- * @brief Floating-point non equivalence comparator based on ULP (Unit in the last place).
- *
- * @note to compare if two floating points `flhs` and `frhs` are not equivalent,
- * use not_float_equal(flhs, frhs), instead of `not_float_equal(flhs-frhs, 0)`.
- * See "Infernal Zero" section of
- * https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
- *
- * @tparam T Type of floating point
- * @tparam max_ulp Maximum tolerable unit in the last place
- * @param flhs First floating point to compare
- * @param frhs Second floating point to compare
- * @return `true` if two floating points differ by greater `ulp`.
- */
-template <typename T, unsigned max_ulp = default_max_ulp>
-bool CUSPATIAL_HOST_DEVICE not_float_equal(T const& flhs, T const& frhs)
-{
-  return !float_equal(flhs, frhs);
-}
-
 }  // namespace detail
 }  // namespace cuspatial

cpp/include/cuspatial/detail/utility/linestring.cuh

@@ -152,7 +152,7 @@ __forceinline__ thrust::optional<segment<T>> __device__ collinear_or_parallel_ov
   auto ac = c - a;
 
   // Parallel
-  if (not_float_equal(det(ab, ac), T{0})) return thrust::nullopt;
+  if (not float_equal(det(ab, ac), T{0})) return thrust::nullopt;
 
   // Must be on the same line, sort the endpoints
   if (b < a) thrust::swap(a, b);
@@ -213,5 +213,21 @@ segment_intersection(segment<T> const& segment1, segment<T> const& segment2)
   return {thrust::nullopt, thrust::nullopt};
 }
 
+/**
+ * @brief Given a point and a segment, returns true if point is on the segment.
+ */
+template <typename T>
+bool __device__ is_point_on_segment(segment<T> const& segment, vec_2d<T> const& c)
+{
+  auto [a, b] = segment;
+  auto ab     = b - a;
+  auto ac     = c - a;
+
+  if (not float_equal(det(ab, ac), T{0})) return false;
+
+  if (b < a) thrust::swap(a, b);
+  return a <= c && c <= b;
+}
+
 }  // namespace detail
 }  // namespace cuspatial

cpp/include/cuspatial/experimental/detail/find/find_point_on_segments.cuh

@@ -0,0 +1,94 @@
+/*
+ * Copyright (c) 2022, 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.
+ */
+
+#pragma once
+
+#include <cuspatial_test/test_util.cuh>
+
+#include <cuspatial/cuda_utils.hpp>
+#include <cuspatial/detail/iterator.hpp>
+#include <cuspatial/detail/utility/linestring.cuh>
+#include <cuspatial/error.hpp>
+
+#include <rmm/cuda_stream_view.hpp>
+#include <rmm/exec_policy.hpp>
+
+#include <thrust/tabulate.h>
+#include <thrust/tuple.h>
+
+namespace cuspatial {
+namespace detail {
+
+/**
+ * @brief Functor to find if the given point is on any of the segments in the same pair
+ */
+template <typename MultiPointRange, typename OffsetsRange, typename SegmentsRange>
+struct find_point_on_segment_functor {
+  MultiPointRange multipoints;
+  OffsetsRange segment_offsets;
+  SegmentsRange segments;
+
+  find_point_on_segment_functor(MultiPointRange multipoints,
+                                OffsetsRange segment_offsets,
+                                SegmentsRange segments)
+    : multipoints(multipoints), segment_offsets(segment_offsets), segments(segments)
+  {
+  }
+
+  template <typename IndexType>
+  uint8_t __device__ operator()(IndexType i)
+  {
+    auto point        = thrust::raw_reference_cast(multipoints.point_begin()[i]);
+    auto geometry_idx = multipoints.geometry_idx_from_point_idx(i);
+
+    for (auto segment_idx = segment_offsets[geometry_idx];
+         segment_idx < segment_offsets[geometry_idx + 1];
+         segment_idx++) {
+      auto const& segment = thrust::raw_reference_cast(segments[segment_idx]);
+      if (is_point_on_segment(segment, point)) return true;
+    }
+    return false;
+  }
+};
+
+/**
+ * @brief Given a multipoint and a set of segments, for each point, if the point is
+ * on any of the segments, set the `mergeable_flag` of the point to `1`.
+ */
+template <typename MultiPointRange,
+          typename OffsetsRange,
+          typename SegmentsRange,
+          typename OutputIt1>
+void find_point_on_segments(MultiPointRange multipoints,
+                            OffsetsRange segment_offsets,
+                            SegmentsRange segments,
+                            OutputIt1 mergeable_flag,
+                            rmm::cuda_stream_view stream)
+{
+  using index_t = typename MultiPointRange::index_t;
+  CUSPATIAL_EXPECTS(multipoints.size() == segments.size(),
+                    "Input should contain the same number of pairs.");
+
+  if (segments.size() == 0) return;
+
+  thrust::tabulate(rmm::exec_policy(stream),
+                   mergeable_flag,
+                   mergeable_flag + multipoints.num_points(),
+                   find_point_on_segment_functor{multipoints, segment_offsets, segments});
+}
+
+}  // namespace detail
+}  // namespace cuspatial

cpp/tests/CMakeLists.txt

@@ -168,4 +168,5 @@ ConfigureTest(OPERATOR_TEST_EXP
     experimental/operators/linestrings_test.cu)
 
 ConfigureTest(FIND_TEST_EXP
+    experimental/find/find_point_on_segments_test.cu
     experimental/find/find_duplicate_points_test.cu)

cpp/tests/experimental/find/find_point_on_segments_test.cu

@@ -0,0 +1,207 @@
+/*
+ * Copyright (c) 2022, 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/vector_equality.hpp>
+#include <cuspatial_test/vector_factories.cuh>
+#include <tests/base_fixture.hpp>
+
+#include <rmm/cuda_stream_view.hpp>
+#include <rmm/device_vector.hpp>
+
+#include <cuspatial/experimental/detail/find/find_point_on_segments.cuh>
+#include <cuspatial/experimental/ranges/range.cuh>
+#include <cuspatial/vec_2d.hpp>
+
+using namespace cuspatial;
+using namespace cuspatial::detail;
+using namespace cuspatial::test;
+
+template <typename T>
+struct FindPointOnSegmentTest : public BaseFixture {
+  rmm::cuda_stream_view stream() { return rmm::cuda_stream_default; }
+};
+
+using TestTypes = ::testing::Types<float, double>;
+TYPED_TEST_CASE(FindPointOnSegmentTest, TestTypes);
+
+TYPED_TEST(FindPointOnSegmentTest, Simple1)
+{
+  using T       = TypeParam;
+  using index_t = std::size_t;
+  using P       = vec_2d<T>;
+  using S       = segment<T>;
+
+  auto multipoints = make_multipoints_array({{P{0.0, 0.0}, P{1.0, 0.0}}});
+
+  auto segment_offsets = make_device_vector<index_t>({0, 1});
+  auto segments        = make_device_vector<segment<T>>({S{P{1.0, 1.0}, P{1.0, -1.0}}});
+
+  rmm::device_vector<uint8_t> flags(multipoints.range().num_points());
+  std::vector<uint8_t> expected_flags{0, 1};
+
+  find_point_on_segments(multipoints.range(),
+                         range(segment_offsets.begin(), segment_offsets.end()),
+                         range(segments.begin(), segments.end()),
+                         flags.begin(),
+                         this->stream());
+
+  CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(flags, expected_flags);
+}
+
+TYPED_TEST(FindPointOnSegmentTest, Simple2)
+{
+  using T       = TypeParam;
+  using index_t = std::size_t;
+  using P       = vec_2d<T>;
+  using S       = segment<T>;
+
+  auto multipoints = make_multipoints_array({{P{0.0, 0.0}, P{1.0, 1.0}}});
+
+  auto segment_offsets = make_device_vector<index_t>({0, 1});
+  auto segments        = make_device_vector<segment<T>>({S{P{2.0, 0.0}, P{0.0, 2.0}}});
+
+  rmm::device_vector<uint8_t> flags(multipoints.range().num_points());
+  std::vector<uint8_t> expected_flags{0, 1};
+
+  find_point_on_segments(multipoints.range(),
+                         range(segment_offsets.begin(), segment_offsets.end()),
+                         range(segments.begin(), segments.end()),
+                         flags.begin(),
+                         this->stream());
+
+  CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(flags, expected_flags);
+}
+
+TYPED_TEST(FindPointOnSegmentTest, Simple3)
+{
+  using T       = TypeParam;
+  using index_t = std::size_t;
+  using P       = vec_2d<T>;
+  using S       = segment<T>;
+
+  auto multipoints = make_multipoints_array({{P{0.0, 0.0}, P{1.0, 1.0}}});
+
+  auto segment_offsets = make_device_vector<index_t>({0, 1});
+  auto segments        = make_device_vector<segment<T>>({S{P{0.0, 1.0}, P{1.0, 1.0}}});
+
+  rmm::device_vector<uint8_t> flags(multipoints.range().num_points());
+  std::vector<uint8_t> expected_flags{0, 1};
+
+  find_point_on_segments(multipoints.range(),
+                         range(segment_offsets.begin(), segment_offsets.end()),
+                         range(segments.begin(), segments.end()),
+                         flags.begin(),
+                         this->stream());
+
+  CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(flags, expected_flags);
+}
+
+TYPED_TEST(FindPointOnSegmentTest, Simple4)
+{
+  using T       = TypeParam;
+  using index_t = std::size_t;
+  using P       = vec_2d<T>;
+  using S       = segment<T>;
+
+  auto multipoints = make_multipoints_array({{P{0.0, 0.0}, P{1.0, 1.0}}});
+
+  auto segment_offsets = make_device_vector<index_t>({0, 1});
+  auto segments        = make_device_vector<segment<T>>({S{P{1.0, 1.0}, P{1.0, 0.0}}});
+
+  rmm::device_vector<uint8_t> flags(multipoints.range().num_points());
+  std::vector<uint8_t> expected_flags{0, 1};
+
+  find_point_on_segments(multipoints.range(),
+                         range(segment_offsets.begin(), segment_offsets.end()),
+                         range(segments.begin(), segments.end()),
+                         flags.begin(),
+                         this->stream());
+
+  CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(flags, expected_flags);
+}
+
+TYPED_TEST(FindPointOnSegmentTest, NoPointOnSegment1)
+{
+  using T       = TypeParam;
+  using index_t = std::size_t;
+  using P       = vec_2d<T>;
+  using S       = segment<T>;
+
+  auto multipoints = make_multipoints_array({{P{0.0, 0.0}, P{1.0, 1.0}}});
+
+  auto segment_offsets = make_device_vector<index_t>({0, 1});
+  auto segments        = make_device_vector<segment<T>>({S{P{0.0, 0.5}, P{1.0, 0.0}}});
+
+  rmm::device_vector<uint8_t> flags(multipoints.range().num_points());
+  std::vector<uint8_t> expected_flags{0, 0};
+
+  find_point_on_segments(multipoints.range(),
+                         range(segment_offsets.begin(), segment_offsets.end()),
+                         range(segments.begin(), segments.end()),
+                         flags.begin(),
+                         this->stream());
+
+  CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(flags, expected_flags);
+}
+
+TYPED_TEST(FindPointOnSegmentTest, NoPointOnSegment2)
+{
+  using T       = TypeParam;
+  using index_t = std::size_t;
+  using P       = vec_2d<T>;
+  using S       = segment<T>;
+
+  auto multipoints = make_multipoints_array({{P{0.0, 0.0}, P{1.0, 1.0}}});
+
+  auto segment_offsets = make_device_vector<index_t>({0, 1});
+  auto segments        = make_device_vector<segment<T>>({S{P{2.0, 2.0}, P{3.0, 3.0}}});
+
+  rmm::device_vector<uint8_t> flags(multipoints.range().num_points());
+  std::vector<uint8_t> expected_flags{0, 0};
+
+  find_point_on_segments(multipoints.range(),
+                         range(segment_offsets.begin(), segment_offsets.end()),
+                         range(segments.begin(), segments.end()),
+                         flags.begin(),
+                         this->stream());
+
+  CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(flags, expected_flags);
+}
+
+TYPED_TEST(FindPointOnSegmentTest, TwoPairs)
+{
+  using T       = TypeParam;
+  using index_t = std::size_t;
+  using P       = vec_2d<T>;
+  using S       = segment<T>;
+
+  auto multipoints = make_multipoints_array({{P{0.0, 0.0}, P{1.0, 1.0}}, {P{2.0, 2.0}}});
+
+  auto segment_offsets = make_device_vector<index_t>({0, 1, 2});
+  auto segments =
+    make_device_vector<segment<T>>({S{P{2.0, 2.0}, P{3.0, 3.0}}, S{P{1.0, 3.0}, P{3.0, 1.0}}});
+
+  rmm::device_vector<uint8_t> flags(multipoints.range().num_points());
+  std::vector<uint8_t> expected_flags{0, 0, 1};
+
+  find_point_on_segments(multipoints.range(),
+                         range(segment_offsets.begin(), segment_offsets.end()),
+                         range(segments.begin(), segments.end()),
+                         flags.begin(),
+                         this->stream());
+
+  CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(flags, expected_flags);
+}