Introduce Segment Intersection Primitive (#778)

Closes #763 , this PR introduces device primitive for segment-segment intersection.

Authors:
  - Michael Wang (https://github.com/isVoid)

Approvers:
  - Mark Harris (https://github.com/harrism)

URL: #778

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

@@ -103,10 +103,15 @@ signmagnitude_to_biased(Bits const& sam)
 /**
  * @brief Floating-point equivalence comparator based on ULP (Unit in the last place).
  *
+ * @note to compare if two floating points `flhs` and `frhs` are equivalent,
+ * use float_equal(flhs, frhs), instead of `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 lhs First floating point to compare
- * @param rhs Second floating point to compare
+ * @param flhs First floating point to compare
+ * @param frhs Second floating point to compare
  * @return `true` if two floating points differ by less or equal to `ulp`.
  */
 template <typename T, unsigned max_ulp = default_max_ulp>
@@ -125,17 +130,21 @@ bool CUSPATIAL_HOST_DEVICE float_equal(T const& flhs, T const& frhs)
 /**
  * @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 lhs First floating point to compare
- * @param rhs Second floating point to compare
+ * @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(FloatingPointBits<T> const& lhs,
-                                           FloatingPointBits<T> const& rhs)
+bool CUSPATIAL_HOST_DEVICE not_float_equal(T const& flhs, T const& frhs)
 {
-  return !float_equal(lhs, rhs);
+  return !float_equal(flhs, frhs);
 }
 
 }  // namespace detail

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

@@ -16,12 +16,18 @@
 
 #pragma once
 
-#include <thrust/tuple.h>
-
+#include <cuspatial/detail/utility/floating_point.cuh>
+#include <cuspatial/experimental/geometry/segment.cuh>
 #include <cuspatial/vec_2d.hpp>
 
+#include <thrust/optional.h>
+#include <thrust/pair.h>
+#include <thrust/swap.h>
+#include <thrust/tuple.h>
+
 namespace cuspatial {
 namespace detail {
+
 /**
  * @internal
  * @brief Get the index that is one-past the end point of linestring at @p linestring_idx
@@ -54,7 +60,7 @@ point_to_segment_distance_squared_nearest_point(vec_2d<T> const& c,
   auto ab        = b - a;
   auto ac        = c - a;
   auto l_squared = dot(ab, ab);
-  if (l_squared == 0) { return thrust::make_tuple(dot(ac, ac), a); }
+  if (float_equal(l_squared, T{0})) { return thrust::make_tuple(dot(ac, ac), a); }
   auto r  = dot(ac, ab);
   auto bc = c - b;
   // If the projection of `c` is outside of segment `ab`, compute point-point distance.
@@ -115,7 +121,7 @@ __forceinline__ T __device__ squared_segment_distance(vec_2d<T> const& a,
   auto cd    = d - c;
   auto denom = det(ab, cd);
 
-  if (denom == 0) {
+  if (float_equal(denom, T{0})) {
     // Segments parallel or collinear
     return segment_distance_no_intersect_or_colinear(a, b, c, d);
   }
@@ -129,5 +135,83 @@ __forceinline__ T __device__ squared_segment_distance(vec_2d<T> const& a,
   return segment_distance_no_intersect_or_colinear(a, b, c, d);
 }
 
+/**
+ * @internal
+ * @brief Given two collinear or parallel segments, return their potential overlapping segment
+ *
+ * @p a, @p b, @p c, @p d refer to end points of segment ab and cd.
+ * @p center is the geometric center of the segments, used to decondition the coordinates.
+ *
+ * @return optional end points of overlapping segment
+ */
+template <typename T>
+__forceinline__ thrust::optional<segment<T>> __device__ collinear_or_parallel_overlapping_segments(
+  vec_2d<T> a, vec_2d<T> b, vec_2d<T> c, vec_2d<T> d, vec_2d<T> center = vec_2d<T>{})
+{
+  auto ab = b - a;
+  auto ac = c - a;
+
+  // Parallel
+  if (not_float_equal(det(ab, ac), T{0})) return thrust::nullopt;
+
+  // Must be on the same line, test if intersect
+  if ((a < c && c < b) || (a < d && d < b)) {
+    // Compute smallest interval between the segments
+    if (b < a) thrust::swap(a, b);
+    if (d < c) thrust::swap(c, d);
+    auto e0 = a > c ? a : c;
+    auto e1 = b < d ? b : d;
+
+    // Decondition the coordinates
+    return segment<T>{e0 + center, e1 + center};
+  }
+
+  return thrust::nullopt;
+}
+
+/**
+ * @internal
+ * @brief Primitive to compute intersections between two segments
+ * Two segments can intersect at a point, overlap at a segment, or be disjoint.
+ *
+ * @return A pair of optional intersecting point and optional overlapping segment
+ */
+template <typename T>
+__forceinline__ thrust::pair<thrust::optional<vec_2d<T>>, thrust::optional<segment<T>>> __device__
+segment_intersection(segment<T> const& segment1, segment<T> const& segment2)
+{
+  auto [a, b] = segment1;
+  auto [c, d] = segment2;
+
+  // Condition the coordinates to avoid large floating point error
+  auto center = midpoint(midpoint(a, b), midpoint(c, d));
+  a -= center;
+  b -= center;
+  c -= center;
+  d -= center;
+
+  auto ab = b - a;
+  auto cd = d - c;
+
+  auto denom = det(ab, cd);
+
+  if (float_equal(denom, T{0})) {
+    // Segments parallel or collinear
+    return {thrust::nullopt, collinear_or_parallel_overlapping_segments(a, b, c, d, center)};
+  }
+
+  auto ac               = c - a;
+  auto r_numer          = det(ac, cd);
+  auto denom_reciprocal = 1 / denom;
+  auto r                = r_numer * denom_reciprocal;
+  auto s                = det(ac, ab) * denom_reciprocal;
+  if (r >= 0 and r <= 1 and s >= 0 and s <= 1) {
+    auto p = a + r * ab;
+    // Decondition the coordinates
+    return {p + center, thrust::nullopt};
+  }
+  return {thrust::nullopt, thrust::nullopt};
+}
+
 }  // namespace detail
 }  // namespace cuspatial

cpp/include/cuspatial/experimental/geometry/segment.cuh

@@ -0,0 +1,30 @@
+/*
+ * 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/vec_2d.hpp>
+
+namespace cuspatial {
+
+template <typename T>
+struct segment {
+  using value_type = T;
+  vec_2d<T> first;
+  vec_2d<T> second;
+};
+
+}  // namespace cuspatial

cpp/include/cuspatial/traits.hpp

@@ -18,7 +18,10 @@
 
 #include <cuspatial/vec_2d.hpp>
 
+#include <thrust/optional.h>
+
 #include <iterator>
+#include <optional>
 #include <type_traits>
 
 namespace cuspatial {
@@ -84,6 +87,13 @@ constexpr bool is_same_floating_point()
          std::conjunction_v<std::is_floating_point<Ts>...>;
 }
 
+template <typename>
+constexpr bool is_optional_impl = false;
+template <typename T>
+constexpr bool is_optional_impl<std::optional<T>> = true;
+template <typename T>
+constexpr bool is_optional = is_optional_impl<std::remove_cv_t<std::remove_reference_t<T>>>;
+
 /**
  * @internal
  * @brief Get the value type of @p Iterator type