Skip to content

Commit

Permalink
Add Multi-Geometry support to point_linestring_distance and build p…
Browse files Browse the repository at this point in the history 
…ython bindings (#660)

This PR adds multi-point and multi-linestring support to `point_linestring_distance`. The c++ API **always** supports multi-variant of the geometry type, since any single-geometry type can be wrapped with a `counting_iteartor` for its underlying geometry offset.

The column API accepts a `std::optional` for the geometry offset inputs. If not provided, they are considered single geometry by default. A `multigeom_dispatcher` is used to generalize the runtime to compile optional information.

This PR also builds python bindings for `point_linestring_distance`. First, the cython bindings is created with a backported `optional` module from cython 3.0 to support the above `std::optional` API. The cython APIs also made use of the dynamic typing property of python to allow geometry offsets to be optional arguments. The python API is the first examplar of a computing API accepting GeoSeries as input. For simplicity, we assume `GeoSeries` only contains single type geometry and there is no mixing of `points` and `multipoints`.

Contributes to #231 
Follow up to #573

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

Approvers:
  - Mark Harris (https://github.com/harrism)
  - H. Thomson Comer (https://github.com/thomcom)

URL: #660
  • Loading branch information
@isVoid
isVoid authored Sep 13, 2022
1 parent e93118a commit 85c3150
Show file tree
Hide file tree
Showing 19 changed files with 1,165 additions and 156 deletions.
51 changes: 43 additions & 8 deletions cpp/include/cuspatial/distance/point_linestring_distance.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,30 +35,65 @@ namespace cuspatial {
* Point: (1, 1)
* Linestring: (0, 0) -> (1, 1) -> (2, 0) -> (3, 0) -> (3, 1)
*
* The input of the abbove example is:
* The input of the above example is:
* multipoint_geometry_offsets: nullopt
* points_xy: {0, 1, 0, 1}
* linestring_offsets: {0, 3, 8}
* multilinestring_geometry_offsets: nullopt
* linestring_part_offsets: {0, 3, 8}
* linestring_xy: {0, 1, 1, 0, 2, 0, 0, 0, 1, 1, 2, 0, 3, 0, 3, 1}
*
* Result: {sqrt(2)/2, 0}
* ```
*
* The following example contains 3 pairs of MultiPoint and MultiLinestring.
* ```
* First pair:
* MultiPoint: (0, 1)
* MultiLinestring: (0, -1) -> (-2, -3), (-4, -5) -> (-5, -6)
*
* Second pair:
* MultiPoint: (2, 3), (4, 5)
* MultiLinestring: (7, 8) -> (8, 9)
*
* Third pair:
* MultiPoint: (6, 7), (8, 9)
* MultiLinestring: (9, 10) -> (10, 11)
* The input of the above example is:
* multipoint_geometry_offsets: {0, 1, 3, 5}
* points_xy: {0, 1, 2, 3, 4, 5, 6, 7, 8, 9}
* multilinestring_geometry_offsets: {0, 2, 3, 5}
* linestring_part_offsets: {0, 2, 4, 6, 8}
* linestring_points_xy: {0, -1, -2, -3, -4, -5, -5, -6, 7, 8, 8, 9, 9, 10, 10 ,11}
*
* Result: {2.0, 4.24264, 1.41421}
* ```
*
* @param multipoint_geometry_offsets Beginning and ending indices to each geometry in the
* multi-point
* @param points_xy Interleaved x, y-coordinates of points
* @param linestring_offsets Beginning and ending indices for each linestring in the `linestring_x`
* and `linestring_y` arrays. Must satisfy `linestring_offsets.size() + 1 == points_xy.size()`.
* @param multilinestring_geometry_offsets Beginning and ending indices to each geometry in the
* multi-linestring
* @param linestring_part_offsets Beginning and ending indices for each linestring in the point
* array. Because the coordinates are interleaved, the actual starting position for the coordinate
* of linestring `i` is `2*linestring_part_offsets[i]`.
* @param linestring_points_xy Interleaved x, y-coordinates of linestring points.
* @param mr Device memory resource used to allocate the returned column.
* @return A column containing the distance between each pair of corresponding points and
* linestrings.
*
* @throws cuspatial::logic_error if the number of points and linestrings do not match.
* @throws cuspatial::logic_error if there is a size mismatch between the x- and y-coordinates of
* the points or linestring points.
* @note Any optional geometry indices, if is `nullopt`, implies the underlying geometry contains
* only one component. Otherwise, it contains multiple components.
*
* @throws cuspatial::logic_error if the number of (multi)points and (multi)linestrings do not
* match.
* @throws cuspatial::logic_error if the any of the point arrays have mismatched types.
*/
std::unique_ptr<cudf::column> pairwise_point_linestring_distance(
std::optional<cudf::device_span<cudf::size_type const>> multipoint_geometry_offsets,
cudf::column_view const& points_xy,
cudf::device_span<cudf::size_type const> linestring_offsets,
std::optional<cudf::device_span<cudf::size_type const>> multilinestring_geometry_offsets,
cudf::device_span<cudf::size_type const> linestring_part_offsets,
cudf::column_view const& linestring_points_xy,
rmm::mr::device_memory_resource* mr = rmm::mr::get_current_device_resource());

Expand Down
143 changes: 104 additions & 39 deletions cpp/include/cuspatial/experimental/detail/point_linestring_distance.cuh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,38 +42,60 @@ namespace detail {

/**
* @internal
* @brief The kernel to compute point to linestring distance
* @brief The kernel to compute multi-point to multi-linestring distance
*
* Each thread computes the distance between a line segment in the linestring and the
* corresponding point in the pair. The shortest distance is computed in the output
* array via an atomic operation.
* corresponding multi-point part in the pair. The shortest distance is computed in the
* output array via an atomic operation.
*
* @tparam Cart2dItA Iterator to 2d cartesian coordinates. Must meet requirements of
* [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
* @tparam Cart2dItB Iterator to 2d cartesian coordinates. Must meet requirements of
* [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
* @tparam OffsetIterator Iterator to linestring offsets. Must meet requirements of
* @tparam OffsetIteratorA Iterator to offsets. Must meet requirements of
* [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
* @tparam OffsetIteratorB Iterator to offsets. Must meet requirements of
* [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
* @tparam OffsetIteratorC Iterator to offsets. Must meet requirements of
* [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
* @tparam OutputIterator Iterator to output distances. Must meet requirements of
* [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible and mutable.
*
* @param[in] point_geometry_offset_first Iterator to the beginning of the range of the multipoint
* parts
* @param[in] point_geometry_offset_last Iterator to the end of the range of the multipoint parts
* @param[in] points_first Iterator to the beginning of the range of the points
* @param[in] points_last Iterator to the end of the range of the points
* @param[in] linestring_offsets_begin Iterator to the beginning of the range of the linestring
* @param[in] linestring_geometry_offset_first Iterator to the beginning of the range of the
* linestring parts
* @param[in] linestring_geometry_offset_last Iterator to the end of the range of the linestring
* parts
* @param[in] linestring_part_offsets_first Iterator to the beginning of the range of the linestring
* offsets
* @param[in] linestring_offsets_end Iterator to the end of the range of the linestring offsets
* @param[in] linestring_points_begin Iterator to the beginning of the range of the linestring
* @param[in] linestring_part_offsets_last Iterator to the beginning of the range of the linestring
* offsets
* @param[in] linestring_points_first Iterator to the beginning of the range of the linestring
* points
* @param[in] linestring_points_end Iterator to the end of the range of the linestring points
* @param[in] linestring_points_last Iterator to the end of the range of the linestring points
* @param[out] distances Iterator to the output range of shortest distances between pairs.
*
* [LinkLRAI]: https://en.cppreference.com/w/cpp/named_req/RandomAccessIterator
* "LegacyRandomAccessIterator"
*/
template <typename Cart2dItA, typename Cart2dItB, typename OffsetIterator, typename OutputIterator>
void __global__ pairwise_point_linestring_distance(Cart2dItA points_first,
OffsetIterator linestring_offsets_first,
OffsetIterator linestring_offsets_last,
template <class Cart2dItA,
class Cart2dItB,
class OffsetIteratorA,
class OffsetIteratorB,
class OffsetIteratorC,
class OutputIterator>
void __global__ pairwise_point_linestring_distance(OffsetIteratorA point_geometry_offset_first,
OffsetIteratorA point_geometry_offset_last,
Cart2dItA points_first,
Cart2dItA points_last,
OffsetIteratorB linestring_geometry_offset_first,
OffsetIteratorB linestring_geometry_offset_last,
OffsetIteratorC linestring_part_offsets_first,
OffsetIteratorC linestring_part_offsets_last,
Cart2dItB linestring_points_first,
Cart2dItB linestring_points_last,
OutputIterator distances)
Expand All @@ -83,34 +105,69 @@ void __global__ pairwise_point_linestring_distance(Cart2dItA points_first,
for (auto idx = threadIdx.x + blockIdx.x * blockDim.x;
idx < std::distance(linestring_points_first, thrust::prev(linestring_points_last));
idx += gridDim.x * blockDim.x) {
auto offsets_iter =
thrust::upper_bound(thrust::seq, linestring_offsets_first, linestring_offsets_last, idx);
// Search from the part offsets array to determine the part idx of current linestring point
auto linestring_part_offsets_iter = thrust::upper_bound(
thrust::seq, linestring_part_offsets_first, linestring_part_offsets_last, idx);

// Pointer to the last point in the linestring, skip iteration.
// Note that the last point for the last linestring is guarded by the grid-stride loop.
if (offsets_iter != linestring_offsets_last and *offsets_iter - 1 == idx) { continue; }

auto point_idx = thrust::distance(linestring_offsets_first, thrust::prev(offsets_iter));
vec_2d<T> const a = linestring_points_first[idx];
vec_2d<T> const b = linestring_points_first[idx + 1];
vec_2d<T> const c = points_first[point_idx];

auto const distance_squared = point_to_segment_distance_squared(c, a, b);

atomicMin(&thrust::raw_reference_cast(*(distances + point_idx)),
static_cast<T>(std::sqrt(distance_squared)));
if (linestring_part_offsets_iter != linestring_part_offsets_last &&
*linestring_part_offsets_iter - 1 == idx) {
continue;
}

auto part_offsets_idx =
thrust::distance(linestring_part_offsets_first, thrust::prev(linestring_part_offsets_iter));

// Search from the linestring geometry offsets array to determine the geometry idx of current
// linestring point
auto geometry_offsets_iter = thrust::upper_bound(thrust::seq,
linestring_geometry_offset_first,
linestring_geometry_offset_last,
part_offsets_idx);
// geometry_idx is also the index to corresponding multipoint in the pair
auto geometry_idx =
thrust::distance(linestring_geometry_offset_first, thrust::prev(geometry_offsets_iter));

// Reduce the minimum distance between different parts of the multi-point.
vec_2d<T> const a = linestring_points_first[idx];
vec_2d<T> const b = linestring_points_first[idx + 1];
T min_distance_squared = std::numeric_limits<T>::max();

for (auto point_idx = point_geometry_offset_first[geometry_idx];
point_idx < point_geometry_offset_first[geometry_idx + 1];
point_idx++) {
vec_2d<T> const c = points_first[point_idx];

// TODO: reduce redundant computation only related to `a`, `b` in this helper.
auto const distance_squared = point_to_segment_distance_squared(c, a, b);
min_distance_squared = std::min(distance_squared, min_distance_squared);
}

atomicMin(&thrust::raw_reference_cast(*(distances + geometry_idx)),
static_cast<T>(std::sqrt(min_distance_squared)));
}
}

} // namespace detail

template <class Cart2dItA, class Cart2dItB, class OffsetIterator, class OutputIt>
void pairwise_point_linestring_distance(Cart2dItA points_first,
Cart2dItA points_last,
OffsetIterator linestring_offsets_first,
Cart2dItB linestring_points_first,
Cart2dItB linestring_points_last,
OutputIt distances_first,
rmm::cuda_stream_view stream)
template <class Cart2dItA,
class Cart2dItB,
class OffsetIteratorA,
class OffsetIteratorB,
class OffsetIteratorC,
class OutputIt>
OutputIt pairwise_point_linestring_distance(OffsetIteratorA point_geometry_offset_first,
OffsetIteratorA point_geometry_offset_last,
Cart2dItA points_first,
Cart2dItA points_last,
OffsetIteratorB linestring_geometry_offset_first,
OffsetIteratorC linestring_part_offsets_first,
OffsetIteratorC linestring_part_offsets_last,
Cart2dItB linestring_points_first,
Cart2dItB linestring_points_last,
OutputIt distances_first,
rmm::cuda_stream_view stream)
{
using T = detail::iterator_vec_base_type<Cart2dItA>;

Expand All @@ -120,11 +177,12 @@ void pairwise_point_linestring_distance(Cart2dItA points_first,
static_assert(detail::is_same<vec_2d<T>,
detail::iterator_value_type<Cart2dItA>,
detail::iterator_value_type<Cart2dItB>>(),
"Inputs must be cuspatial::vec_2d<T>");
"Inputs must be cuspatial::vec_2d");

auto const num_pairs = thrust::distance(points_first, points_last);
auto const num_pairs =
thrust::distance(point_geometry_offset_first, point_geometry_offset_last) - 1;

if (num_pairs == 0) { return; }
if (num_pairs == 0) { return distances_first; }

auto const num_linestring_points =
thrust::distance(linestring_points_first, linestring_points_last);
Expand All @@ -137,14 +195,21 @@ void pairwise_point_linestring_distance(Cart2dItA points_first,
(num_linestring_points + threads_per_block - 1) / threads_per_block;

detail::pairwise_point_linestring_distance<<<num_blocks, threads_per_block, 0, stream.value()>>>(
point_geometry_offset_first,
point_geometry_offset_last,
points_first,
linestring_offsets_first,
linestring_offsets_first + num_pairs + 1,
points_last,
linestring_geometry_offset_first,
linestring_geometry_offset_first + num_pairs + 1,
linestring_part_offsets_first,
linestring_part_offsets_last,
linestring_points_first,
linestring_points_first + num_linestring_points,
linestring_points_last,
distances_first);

CUSPATIAL_CUDA_TRY(cudaGetLastError());

return distances_first + num_pairs;
}

} // namespace cuspatial
56 changes: 38 additions & 18 deletions cpp/include/cuspatial/experimental/point_linestring_distance.cuh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,20 +30,30 @@ namespace cuspatial {
* the requirements of [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
* @tparam Cart2dItB iterator type for point array of the linestring element of each pair. Must meet
* the requirements of [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
* @tparam OffsetIterator iterator type for offset array. Must meet the requirements of
* [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
* @tparam OffsetIteratorA iterator type for `point_geometry_offset` array. Must meet the
* requirements of [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
* @tparam OffsetIteratorB iterator type for `linestring_geometry_offset` array. Must meet the
* requirements of [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
* @tparam OffsetIteratorC iterator type for `linestring_part_offset` array. Must meet the
* requirements of [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
* @tparam OutputIt iterator type for output array. Must meet the requirements of
* [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
*
* @param points_first beginning of the range of points making up the first element of each
* pair
* @param points_last end of the range of the points making up the first element of each pair
* @param linestring_offsets_first beginning of the range of the offsets to the linestring element
* of each pair
* @param linestring_points_first beginning of the range of points of the linestring element of each
* @param point_geometry_offset_first beginning of the range of multipoint geometries of each
* pair
* @param linestring_points_last end of the range of points of the linestring element of each pair
* @param distances_first beginning the output range of distances
* @param point_geometry_offset_last end of the range of multipoint geometries of each pair
* @param points_first beginning of the range of point values
* @param points_last end of the range of the point values
* @param linestring_geometry_offset_first beginning of the range of offsets to the multilinestring
* geometry of each pair, the end range is implied by linestring_geometry_offset_first +
* std::distance(`point_geometry_offset_first`, `point_geometry_offset_last`)
* @param linestring_offsets_first beginning of the range of offsets to the starting point
* of each linestring
* @param linestring_offsets_last end of the range of offsets to the starting point
* of each linestring
* @param linestring_points_first beginning of the range of linestring points
* @param linestring_points_last end of the range of linestring points
* @param distances_first beginning of the output range of distances
* @param stream The CUDA stream to use for device memory operations and kernel launches.
*
* @pre all input iterators for coordinates must have `cuspatial::vec_2d` type.
Expand All @@ -53,14 +63,24 @@ namespace cuspatial {
* [LinkLRAI]: https://en.cppreference.com/w/cpp/named_req/RandomAccessIterator
* "LegacyRandomAccessIterator"
*/
template <class Cart2dItA, class Cart2dItB, class OffsetIterator, class OutputIt>
void pairwise_point_linestring_distance(Cart2dItA points_first,
Cart2dItA points_last,
OffsetIterator linestring_offsets_first,
Cart2dItB linestring_points_first,
Cart2dItB linestring_points_last,
OutputIt distances_first,
rmm::cuda_stream_view stream = rmm::cuda_stream_default);
template <class Cart2dItA,
class Cart2dItB,
class OffsetIteratorA,
class OffsetIteratorB,
class OffsetIteratorC,
class OutputIt>
OutputIt pairwise_point_linestring_distance(
OffsetIteratorA point_geometry_offset_first,
OffsetIteratorA point_geometry_offset_last,
Cart2dItA points_first,
Cart2dItA points_last,
OffsetIteratorB linestring_geometry_offset_first,
OffsetIteratorC linestring_part_offsets_first,
OffsetIteratorC linestring_part_offsets_last,
Cart2dItB linestring_points_first,
Cart2dItB linestring_points_last,
OutputIt distances_first,
rmm::cuda_stream_view stream = rmm::cuda_stream_default);

} // namespace cuspatial

Expand Down
Loading

0 comments on commit 85c3150

Please sign in to comment.