pwrliang commented on code in PR #310: URL: https://github.com/apache/sedona-db/pull/310#discussion_r2591534237
########## c/sedona-libgpuspatial/libgpuspatial/test/related_test.cu: ########## @@ -0,0 +1,1287 @@ +// Licensed to the Apache Software Foundation (ASF) under one +// or more contributor license agreements. See the NOTICE file +// distributed with this work for additional information +// regarding copyright ownership. The ASF licenses this file +// to you 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 "array_stream.hpp" +#include "gpuspatial/loader/parallel_wkb_loader.h" +#include "gpuspatial/relate/relate.cuh" +#include "gpuspatial/utils/pinned_vector.h" + +#include "test_common.hpp" + +#include <rmm/cuda_stream_view.hpp> + +#include <geos/geom/Geometry.h> +#include <geos/io/WKTReader.h> +#include <geos/operation/relateng/RelateGeometry.h> +#include <geos/operation/relateng/RelateMatrixPredicate.h> +#include <geos/operation/relateng/RelateNG.h> +#include <geos/operation/relateng/RelatePredicate.h> +#include <gtest/gtest.h> + +using namespace geos::geom; +using namespace geos::operation::relateng; +using geos::io::WKTReader; + +// Test cases are from +// https://github.com/libgeos/geos/blob/2d2802d7f7acd7919599b94f3d1530e8cd987aee/tests/unit/operation/relateng/RelateNGTest.cpp + +namespace gpuspatial { +using point_t = Point<double, 2>; +using index_t = uint32_t; +using box_t = Box<Point<float, 2>>; +using loader_t = ParallelWkbLoader<point_t, index_t>; + +template <typename POINT_T, typename INDEX_T> +struct Context { + PinnedVector<POINT_T> points; + PinnedVector<INDEX_T> prefix_sum1; + PinnedVector<INDEX_T> prefix_sum2; + PinnedVector<INDEX_T> prefix_sum3; + PinnedVector<box_t> mbrs; +}; + +template <typename POINT_T> +void ParseWKTPoint(const char* wkt, POINT_T& point) { + nanoarrow::UniqueArrayStream stream; + ArrayStreamFromWKT({{wkt}}, GEOARROW_TYPE_WKB, stream.get()); + nanoarrow::UniqueArray array; + ArrowError error; + ArrowErrorSet(&error, ""); + + ASSERT_EQ(ArrowArrayStreamGetNext(stream.get(), array.get(), &error), NANOARROW_OK); + loader_t loader; + auto cuda_stream = rmm::cuda_stream_default; + + loader.Init(); + loader.Parse(cuda_stream, array.get(), 0, array->length); + auto device_geometries = loader.Finish(cuda_stream); + auto h_vec = TestUtils::ToVector(cuda_stream, device_geometries.get_points()); + cuda_stream.synchronize(); + point = h_vec[0]; +} + +template <typename POINT_T, typename INDEX_T> +void ParseWKTMultiPoint(Context<POINT_T, INDEX_T>& ctx, const char* wkt, + MultiPoint<POINT_T>& multi_point) { + nanoarrow::UniqueArrayStream stream; + ArrayStreamFromWKT({{wkt}}, GEOARROW_TYPE_WKB, stream.get()); + nanoarrow::UniqueArray array; + ArrowError error; + ArrowErrorSet(&error, ""); + + ASSERT_EQ(ArrowArrayStreamGetNext(stream.get(), array.get(), &error), NANOARROW_OK); + loader_t loader; + auto cuda_stream = rmm::cuda_stream_default; + + loader.Init(); + loader.Parse(cuda_stream, array.get(), 0, array->length); + auto device_geometries = loader.Finish(cuda_stream); + + ctx.prefix_sum1 = TestUtils::ToVector( + cuda_stream, device_geometries.get_offsets().multi_point_offsets.ps_num_points); + ctx.points = TestUtils::ToVector(cuda_stream, device_geometries.get_points()); + ctx.mbrs = TestUtils::ToVector(cuda_stream, device_geometries.get_mbrs()); + cuda_stream.synchronize(); + MultiPointArrayView multi_array_view( + ArrayView<INDEX_T>(ctx.prefix_sum1.data(), ctx.prefix_sum1.size()), + ArrayView<POINT_T>(ctx.points.data(), ctx.points.size()), + ArrayView<box_t>(ctx.mbrs.data(), ctx.mbrs.size())); + multi_point = multi_array_view[0]; +} + +template <typename POINT_T, typename INDEX_T> +void ParseWKTLineString(Context<POINT_T, INDEX_T>& ctx, const char* wkt, + LineString<POINT_T>& ls) { + nanoarrow::UniqueArrayStream stream; + ArrayStreamFromWKT({{wkt}}, GEOARROW_TYPE_WKB, stream.get()); + nanoarrow::UniqueArray array; + ArrowError error; + ArrowErrorSet(&error, ""); + + ASSERT_EQ(ArrowArrayStreamGetNext(stream.get(), array.get(), &error), NANOARROW_OK); + loader_t loader; + auto cuda_stream = rmm::cuda_stream_default; + + loader.Init(); + loader.Parse(cuda_stream, array.get(), 0, array->length); + auto device_geometries = loader.Finish(cuda_stream); + ctx.prefix_sum1 = TestUtils::ToVector( + cuda_stream, device_geometries.get_offsets().line_string_offsets.ps_num_points); + ctx.points = TestUtils::ToVector(cuda_stream, device_geometries.get_points()); + ctx.mbrs = TestUtils::ToVector(cuda_stream, device_geometries.get_mbrs()); + cuda_stream.synchronize(); + LineStringArrayView<POINT_T, INDEX_T> ls_array_view( + ArrayView<INDEX_T>(ctx.prefix_sum1.data(), ctx.prefix_sum1.size()), + ArrayView<POINT_T>(ctx.points.data(), ctx.points.size()), + ArrayView<box_t>(ctx.mbrs.data(), ctx.mbrs.size())); + ls = ls_array_view[0]; +} + +template <typename POINT_T, typename INDEX_T> +void ParseWKTMultiLineString(Context<POINT_T, INDEX_T>& ctx, const char* wkt, + MultiLineString<POINT_T, INDEX_T>& m_ls) { + nanoarrow::UniqueArrayStream stream; + ArrayStreamFromWKT({{wkt}}, GEOARROW_TYPE_WKB, stream.get()); + nanoarrow::UniqueArray array; + ArrowError error; + ArrowErrorSet(&error, ""); + + ASSERT_EQ(ArrowArrayStreamGetNext(stream.get(), array.get(), &error), NANOARROW_OK); + loader_t loader; + auto cuda_stream = rmm::cuda_stream_default; + + loader.Init(); + loader.Parse(cuda_stream, array.get(), 0, array->length); + auto device_geometries = loader.Finish(cuda_stream); + ctx.prefix_sum1 = TestUtils::ToVector( + cuda_stream, + device_geometries.get_offsets().multi_line_string_offsets.ps_num_parts); + ctx.prefix_sum2 = TestUtils::ToVector( + cuda_stream, + device_geometries.get_offsets().multi_line_string_offsets.ps_num_points); + ctx.points = TestUtils::ToVector(cuda_stream, device_geometries.get_points()); + ctx.mbrs = TestUtils::ToVector(cuda_stream, device_geometries.get_mbrs()); + cuda_stream.synchronize(); + MultiLineStringArrayView<POINT_T, INDEX_T> m_ls_array_view( + ArrayView<INDEX_T>(ctx.prefix_sum1.data(), ctx.prefix_sum1.size()), + ArrayView<INDEX_T>(ctx.prefix_sum2.data(), ctx.prefix_sum2.size()), + ArrayView<POINT_T>(ctx.points.data(), ctx.points.size()), + ArrayView<box_t>(ctx.mbrs.data(), ctx.mbrs.size())); + m_ls = m_ls_array_view[0]; +} + +template <typename POINT_T, typename INDEX_T> +void ParseWKTPolygon(Context<POINT_T, INDEX_T>& ctx, const char* wkt, + Polygon<POINT_T, INDEX_T>& poly) { + nanoarrow::UniqueArrayStream stream; + ArrayStreamFromWKT({{wkt}}, GEOARROW_TYPE_WKB, stream.get()); + nanoarrow::UniqueArray array; + ArrowError error; + ArrowErrorSet(&error, ""); + + ASSERT_EQ(ArrowArrayStreamGetNext(stream.get(), array.get(), &error), NANOARROW_OK); + loader_t loader; + auto cuda_stream = rmm::cuda_stream_default; + + loader.Init(); + loader.Parse(cuda_stream, array.get(), 0, array->length); + auto device_geometries = loader.Finish(cuda_stream); + ctx.prefix_sum1 = TestUtils::ToVector( + cuda_stream, device_geometries.get_offsets().polygon_offsets.ps_num_rings); + ctx.prefix_sum2 = TestUtils::ToVector( + cuda_stream, device_geometries.get_offsets().polygon_offsets.ps_num_points); + ctx.points = TestUtils::ToVector(cuda_stream, device_geometries.get_points()); + ctx.mbrs = TestUtils::ToVector(cuda_stream, device_geometries.get_mbrs()); + cuda_stream.synchronize(); + PolygonArrayView<POINT_T, INDEX_T> poly_array_view( + ArrayView<INDEX_T>(ctx.prefix_sum1.data(), ctx.prefix_sum1.size()), + ArrayView<INDEX_T>(ctx.prefix_sum2.data(), ctx.prefix_sum2.size()), + ArrayView<POINT_T>(ctx.points.data(), ctx.points.size()), + ArrayView<box_t>(ctx.mbrs.data(), ctx.mbrs.size())); + poly = poly_array_view[0]; +} + +template <typename POINT_T, typename INDEX_T> +void ParseWKTMultiPolygon(Context<POINT_T, INDEX_T>& ctx, const char* wkt, + MultiPolygon<POINT_T, INDEX_T>& poly) { + nanoarrow::UniqueArrayStream stream; + ArrayStreamFromWKT({{wkt}}, GEOARROW_TYPE_WKB, stream.get()); + nanoarrow::UniqueArray array; + ArrowError error; + ArrowErrorSet(&error, ""); + + ASSERT_EQ(ArrowArrayStreamGetNext(stream.get(), array.get(), &error), NANOARROW_OK); + loader_t loader; + auto cuda_stream = rmm::cuda_stream_default; + + loader.Init(); + loader.Parse(cuda_stream, array.get(), 0, array->length); + auto device_geometries = loader.Finish(cuda_stream); + ctx.prefix_sum1 = TestUtils::ToVector( + cuda_stream, device_geometries.get_offsets().multi_polygon_offsets.ps_num_parts); + ctx.prefix_sum2 = TestUtils::ToVector( + cuda_stream, device_geometries.get_offsets().multi_polygon_offsets.ps_num_rings); + ctx.prefix_sum3 = TestUtils::ToVector( + cuda_stream, device_geometries.get_offsets().multi_polygon_offsets.ps_num_points); + ctx.points = TestUtils::ToVector(cuda_stream, device_geometries.get_points()); + ctx.mbrs = TestUtils::ToVector(cuda_stream, device_geometries.get_mbrs()); + cuda_stream.synchronize(); + MultiPolygonArrayView<POINT_T, INDEX_T> poly_array_view( + ArrayView<INDEX_T>(ctx.prefix_sum1.data(), ctx.prefix_sum1.size()), + ArrayView<INDEX_T>(ctx.prefix_sum2.data(), ctx.prefix_sum2.size()), + ArrayView<INDEX_T>(ctx.prefix_sum3.data(), ctx.prefix_sum3.size()), + ArrayView<POINT_T>(ctx.points.data(), ctx.points.size()), + ArrayView<box_t>(ctx.mbrs.data(), ctx.mbrs.size())); + poly = poly_array_view[0]; +} + +template <typename GEOMETRY1_T, typename GEOMETRY2_T> +void TestRelate(const char* wkt1, const char* wkt2, const GEOMETRY1_T& g1, + const GEOMETRY2_T& g2) { + WKTReader r; + auto a = r.read(wkt1); + auto b = r.read(wkt2); + + RelateMatrixPredicate pred; + RelateNG::relate(a.get(), b.get(), pred); + std::string actualVal = pred.getIM()->toString(); + + int val = relate(g1, g2); + char res[10]; + IM__ToString(val, res); + ASSERT_STREQ(actualVal.c_str(), res); +} + +TEST(RelateTest, PointPointDisjoint) { + point_t p1, p2; + + std::string wkt1 = "POINT (0 0)"; + std::string wkt2 = "POINT (1 1)"; + ParseWKTPoint(wkt1.c_str(), p1); + ParseWKTPoint(wkt2.c_str(), p2); + TestRelate(wkt1.c_str(), wkt2.c_str(), p1, p2); +} + +TEST(RelateTest, MultiPointMultiPointContained) { + MultiPoint<point_t> m1, m2; + std::string wkt1 = "MULTIPOINT (0 0, 1 1, 2 2)"; + std::string wkt2 = "MULTIPOINT (1 1, 2 2)"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTMultiPoint(ctx1, wkt1.c_str(), m1); + ParseWKTMultiPoint(ctx2, wkt2.c_str(), m2); + + TestRelate(wkt1.c_str(), wkt2.c_str(), m1, m2); +} + +TEST(RelateTest, MultiPointMultiPointEqual) { + using point_t = Point<double, 2>; + MultiPoint<point_t> m1, m2; + std::string wkt1 = "MULTIPOINT (0 0, 1 1, 2 2)"; + std::string wkt2 = "MULTIPOINT (0 0, 1 1, 2 2)"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTMultiPoint(ctx1, wkt1.c_str(), m1); + ParseWKTMultiPoint(ctx2, wkt2.c_str(), m2); + + TestRelate(wkt1.c_str(), wkt2.c_str(), m1, m2); +} + +TEST(RelateTest, MultiPointMultiPointValidateRelatePP_13) { + MultiPoint<point_t> m1, m2; + std::string wkt1 = "MULTIPOINT ((80 70), (140 120), (20 20), (200 170))"; + std::string wkt2 = "MULTIPOINT ((80 70), (140 120), (80 170), (200 80))"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTMultiPoint(ctx1, wkt1.c_str(), m1); + ParseWKTMultiPoint(ctx2, wkt2.c_str(), m2); + + TestRelate(wkt1.c_str(), wkt2.c_str(), m1, m2); +} + +TEST(RelateTest, LineStringMultiPointContains) { + LineString<point_t> ls1; + MultiPoint<point_t> m2; + std::string wkt1 = "LINESTRING (0 0, 1 1, 2 2)"; + std::string wkt2 = "MULTIPOINT (0 0, 1 1, 2 2)"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTLineString(ctx1, wkt1.c_str(), ls1); + ParseWKTMultiPoint(ctx2, wkt2.c_str(), m2); + TestRelate(wkt1.c_str(), wkt2.c_str(), ls1, m2); +} + +TEST(RelateTest, LineStringMultiPointOverlaps) { + LineString<point_t> ls1; + MultiPoint<point_t> m2; + std::string wkt1 = "LINESTRING (0 0, 1 1)"; + std::string wkt2 = "MULTIPOINT (0 0, 1 1, 2 2)"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTLineString(ctx1, wkt1.c_str(), ls1); + ParseWKTMultiPoint(ctx2, wkt2.c_str(), m2); + TestRelate(wkt1.c_str(), wkt2.c_str(), ls1, m2); +} + +TEST(RelateTest, ZeroLengthLinePoint) { + LineString<point_t> ls1; + point_t p2; + std::string wkt1 = "LINESTRING (0 0, 0 0)"; + std::string wkt2 = "POINT (0 0)"; + Context<point_t, index_t> ctx1; + + ParseWKTLineString(ctx1, wkt1.c_str(), ls1); + ParseWKTPoint(wkt2.c_str(), p2); + TestRelate(wkt1.c_str(), wkt2.c_str(), ls1, p2); +} + +TEST(RelateTest, ZeroLengthLineLine) { + LineString<point_t> ls1; + LineString<point_t> ls2; + std::string wkt1 = "LINESTRING (10 10, 10 10, 10 10)"; + std::string wkt2 = "LINESTRING (10 10, 10 10)"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTLineString(ctx1, wkt1.c_str(), ls1); + ParseWKTLineString(ctx2, wkt2.c_str(), ls2); + TestRelate(wkt1.c_str(), wkt2.c_str(), ls1, ls2); +} + +TEST(RelateTest, NonZeroLengthLinePoint) { + LineString<point_t> ls1; + point_t p2; + std::string wkt1 = "LINESTRING (0 0, 0 0, 9 9)"; + std::string wkt2 = "POINT (1 1)"; + Context<point_t, index_t> ctx1; + + ParseWKTLineString(ctx1, wkt1.c_str(), ls1); + ParseWKTPoint(wkt2.c_str(), p2); + TestRelate(wkt1.c_str(), wkt2.c_str(), ls1, p2); +} + +TEST(RelateTest, LinePointIntAndExt) { + MultiPoint<point_t> m1; + LineString<point_t> ls2; + std::string wkt1 = "MULTIPOINT ((60 60), (100 100))"; + std::string wkt2 = "LINESTRING (40 40, 80 80)"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTMultiPoint(ctx1, wkt1.c_str(), m1); + ParseWKTLineString(ctx2, wkt2.c_str(), ls2); + TestRelate(wkt1.c_str(), wkt2.c_str(), m1, ls2); +} + +TEST(RelateTest, LinesCrossProper) { + LineString<point_t> ls1; + LineString<point_t> ls2; + std::string wkt1 = "LINESTRING (0 0, 9 9)"; + std::string wkt2 = "LINESTRING (0 9, 9 0)"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTLineString(ctx1, wkt1.c_str(), ls1); + ParseWKTLineString(ctx2, wkt2.c_str(), ls2); + TestRelate(wkt1.c_str(), wkt2.c_str(), ls1, ls2); +} + +TEST(RelateTest, LinesOverlap) { + LineString<point_t> ls1; + LineString<point_t> ls2; + std::string wkt1 = "LINESTRING (0 0, 5 5)"; + std::string wkt2 = "LINESTRING (3 3, 9 9)"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTLineString(ctx1, wkt1.c_str(), ls1); + ParseWKTLineString(ctx2, wkt2.c_str(), ls2); + TestRelate(wkt1.c_str(), wkt2.c_str(), ls1, ls2); +} + +TEST(RelateTest, LinesCrossVertex) { + LineString<point_t> ls1; + LineString<point_t> ls2; + std::string wkt1 = "LINESTRING (0 0, 8 8)"; + std::string wkt2 = "LINESTRING (0 8, 4 4, 8 0)"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTLineString(ctx1, wkt1.c_str(), ls1); + ParseWKTLineString(ctx2, wkt2.c_str(), ls2); + TestRelate(wkt1.c_str(), wkt2.c_str(), ls1, ls2); +} + +TEST(RelateTest, LinesTouchVertex) { + LineString<point_t> ls1; + LineString<point_t> ls2; + std::string wkt1 = "LINESTRING (0 0, 8 0)"; + std::string wkt2 = "LINESTRING (0 8, 4 0, 8 8)"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTLineString(ctx1, wkt1.c_str(), ls1); + ParseWKTLineString(ctx2, wkt2.c_str(), ls2); + TestRelate(wkt1.c_str(), wkt2.c_str(), ls1, ls2); +} + +TEST(RelateTest, LinesDisjointByEnvelope) { + LineString<point_t> ls1; + LineString<point_t> ls2; + std::string wkt1 = "LINESTRING (0 0, 9 9)"; + std::string wkt2 = "LINESTRING (10 19, 19 10)"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTLineString(ctx1, wkt1.c_str(), ls1); + ParseWKTLineString(ctx2, wkt2.c_str(), ls2); + TestRelate(wkt1.c_str(), wkt2.c_str(), ls1, ls2); +} + +TEST(RelateTest, LinesDisjoint) { + LineString<point_t> ls1; + LineString<point_t> ls2; + std::string wkt1 = "LINESTRING (0 0, 9 9)"; + std::string wkt2 = "LINESTRING (4 2, 8 6)"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTLineString(ctx1, wkt1.c_str(), ls1); + ParseWKTLineString(ctx2, wkt2.c_str(), ls2); + TestRelate(wkt1.c_str(), wkt2.c_str(), ls1, ls2); +} + +// FIXME: wkt1 is a closed polyline, which has no boundary according to JTS's +// Mod2BoundaryNodeRule We have to implement a similar rule in gpuspatial to handle this +// case correctly TEST(RelateTest, LinesClosedEmpty) { +// MultiLineString<point_t, index_t> m_ls1; +// LineString<point_t> ls2; +// std::string wkt1 = "MULTILINESTRING ((0 0, 0 1), (0 1, 1 1, 1 0, 0 0))"; +// std::string wkt2 = "LINESTRING EMPTY"; +// Context<point_t, index_t> ctx1, ctx2; +// +// ParseWKTMultiLineString(ctx1, wkt1.c_str(), m_ls1); +// ParseWKTLineString(ctx2, wkt2.c_str(), ls2); +// TestRelate(wkt1.c_str(), wkt2.c_str(), m_ls1, ls2); +// } Review Comment: https://github.com/apache/sedona-db/issues/415 ########## c/sedona-libgpuspatial/libgpuspatial/test/related_test.cu: ########## @@ -0,0 +1,1287 @@ +// Licensed to the Apache Software Foundation (ASF) under one +// or more contributor license agreements. See the NOTICE file +// distributed with this work for additional information +// regarding copyright ownership. The ASF licenses this file +// to you 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 "array_stream.hpp" +#include "gpuspatial/loader/parallel_wkb_loader.h" +#include "gpuspatial/relate/relate.cuh" +#include "gpuspatial/utils/pinned_vector.h" + +#include "test_common.hpp" + +#include <rmm/cuda_stream_view.hpp> + +#include <geos/geom/Geometry.h> +#include <geos/io/WKTReader.h> +#include <geos/operation/relateng/RelateGeometry.h> +#include <geos/operation/relateng/RelateMatrixPredicate.h> +#include <geos/operation/relateng/RelateNG.h> +#include <geos/operation/relateng/RelatePredicate.h> +#include <gtest/gtest.h> + +using namespace geos::geom; +using namespace geos::operation::relateng; +using geos::io::WKTReader; + +// Test cases are from +// https://github.com/libgeos/geos/blob/2d2802d7f7acd7919599b94f3d1530e8cd987aee/tests/unit/operation/relateng/RelateNGTest.cpp + +namespace gpuspatial { +using point_t = Point<double, 2>; +using index_t = uint32_t; +using box_t = Box<Point<float, 2>>; +using loader_t = ParallelWkbLoader<point_t, index_t>; + +template <typename POINT_T, typename INDEX_T> +struct Context { + PinnedVector<POINT_T> points; + PinnedVector<INDEX_T> prefix_sum1; + PinnedVector<INDEX_T> prefix_sum2; + PinnedVector<INDEX_T> prefix_sum3; + PinnedVector<box_t> mbrs; +}; + +template <typename POINT_T> +void ParseWKTPoint(const char* wkt, POINT_T& point) { + nanoarrow::UniqueArrayStream stream; + ArrayStreamFromWKT({{wkt}}, GEOARROW_TYPE_WKB, stream.get()); + nanoarrow::UniqueArray array; + ArrowError error; + ArrowErrorSet(&error, ""); + + ASSERT_EQ(ArrowArrayStreamGetNext(stream.get(), array.get(), &error), NANOARROW_OK); + loader_t loader; + auto cuda_stream = rmm::cuda_stream_default; + + loader.Init(); + loader.Parse(cuda_stream, array.get(), 0, array->length); + auto device_geometries = loader.Finish(cuda_stream); + auto h_vec = TestUtils::ToVector(cuda_stream, device_geometries.get_points()); + cuda_stream.synchronize(); + point = h_vec[0]; +} + +template <typename POINT_T, typename INDEX_T> +void ParseWKTMultiPoint(Context<POINT_T, INDEX_T>& ctx, const char* wkt, + MultiPoint<POINT_T>& multi_point) { + nanoarrow::UniqueArrayStream stream; + ArrayStreamFromWKT({{wkt}}, GEOARROW_TYPE_WKB, stream.get()); + nanoarrow::UniqueArray array; + ArrowError error; + ArrowErrorSet(&error, ""); + + ASSERT_EQ(ArrowArrayStreamGetNext(stream.get(), array.get(), &error), NANOARROW_OK); + loader_t loader; + auto cuda_stream = rmm::cuda_stream_default; + + loader.Init(); + loader.Parse(cuda_stream, array.get(), 0, array->length); + auto device_geometries = loader.Finish(cuda_stream); + + ctx.prefix_sum1 = TestUtils::ToVector( + cuda_stream, device_geometries.get_offsets().multi_point_offsets.ps_num_points); + ctx.points = TestUtils::ToVector(cuda_stream, device_geometries.get_points()); + ctx.mbrs = TestUtils::ToVector(cuda_stream, device_geometries.get_mbrs()); + cuda_stream.synchronize(); + MultiPointArrayView multi_array_view( + ArrayView<INDEX_T>(ctx.prefix_sum1.data(), ctx.prefix_sum1.size()), + ArrayView<POINT_T>(ctx.points.data(), ctx.points.size()), + ArrayView<box_t>(ctx.mbrs.data(), ctx.mbrs.size())); + multi_point = multi_array_view[0]; +} + +template <typename POINT_T, typename INDEX_T> +void ParseWKTLineString(Context<POINT_T, INDEX_T>& ctx, const char* wkt, + LineString<POINT_T>& ls) { + nanoarrow::UniqueArrayStream stream; + ArrayStreamFromWKT({{wkt}}, GEOARROW_TYPE_WKB, stream.get()); + nanoarrow::UniqueArray array; + ArrowError error; + ArrowErrorSet(&error, ""); + + ASSERT_EQ(ArrowArrayStreamGetNext(stream.get(), array.get(), &error), NANOARROW_OK); + loader_t loader; + auto cuda_stream = rmm::cuda_stream_default; + + loader.Init(); + loader.Parse(cuda_stream, array.get(), 0, array->length); + auto device_geometries = loader.Finish(cuda_stream); + ctx.prefix_sum1 = TestUtils::ToVector( + cuda_stream, device_geometries.get_offsets().line_string_offsets.ps_num_points); + ctx.points = TestUtils::ToVector(cuda_stream, device_geometries.get_points()); + ctx.mbrs = TestUtils::ToVector(cuda_stream, device_geometries.get_mbrs()); + cuda_stream.synchronize(); + LineStringArrayView<POINT_T, INDEX_T> ls_array_view( + ArrayView<INDEX_T>(ctx.prefix_sum1.data(), ctx.prefix_sum1.size()), + ArrayView<POINT_T>(ctx.points.data(), ctx.points.size()), + ArrayView<box_t>(ctx.mbrs.data(), ctx.mbrs.size())); + ls = ls_array_view[0]; +} + +template <typename POINT_T, typename INDEX_T> +void ParseWKTMultiLineString(Context<POINT_T, INDEX_T>& ctx, const char* wkt, + MultiLineString<POINT_T, INDEX_T>& m_ls) { + nanoarrow::UniqueArrayStream stream; + ArrayStreamFromWKT({{wkt}}, GEOARROW_TYPE_WKB, stream.get()); + nanoarrow::UniqueArray array; + ArrowError error; + ArrowErrorSet(&error, ""); + + ASSERT_EQ(ArrowArrayStreamGetNext(stream.get(), array.get(), &error), NANOARROW_OK); + loader_t loader; + auto cuda_stream = rmm::cuda_stream_default; + + loader.Init(); + loader.Parse(cuda_stream, array.get(), 0, array->length); + auto device_geometries = loader.Finish(cuda_stream); + ctx.prefix_sum1 = TestUtils::ToVector( + cuda_stream, + device_geometries.get_offsets().multi_line_string_offsets.ps_num_parts); + ctx.prefix_sum2 = TestUtils::ToVector( + cuda_stream, + device_geometries.get_offsets().multi_line_string_offsets.ps_num_points); + ctx.points = TestUtils::ToVector(cuda_stream, device_geometries.get_points()); + ctx.mbrs = TestUtils::ToVector(cuda_stream, device_geometries.get_mbrs()); + cuda_stream.synchronize(); + MultiLineStringArrayView<POINT_T, INDEX_T> m_ls_array_view( + ArrayView<INDEX_T>(ctx.prefix_sum1.data(), ctx.prefix_sum1.size()), + ArrayView<INDEX_T>(ctx.prefix_sum2.data(), ctx.prefix_sum2.size()), + ArrayView<POINT_T>(ctx.points.data(), ctx.points.size()), + ArrayView<box_t>(ctx.mbrs.data(), ctx.mbrs.size())); + m_ls = m_ls_array_view[0]; +} + +template <typename POINT_T, typename INDEX_T> +void ParseWKTPolygon(Context<POINT_T, INDEX_T>& ctx, const char* wkt, + Polygon<POINT_T, INDEX_T>& poly) { + nanoarrow::UniqueArrayStream stream; + ArrayStreamFromWKT({{wkt}}, GEOARROW_TYPE_WKB, stream.get()); + nanoarrow::UniqueArray array; + ArrowError error; + ArrowErrorSet(&error, ""); + + ASSERT_EQ(ArrowArrayStreamGetNext(stream.get(), array.get(), &error), NANOARROW_OK); + loader_t loader; + auto cuda_stream = rmm::cuda_stream_default; + + loader.Init(); + loader.Parse(cuda_stream, array.get(), 0, array->length); + auto device_geometries = loader.Finish(cuda_stream); + ctx.prefix_sum1 = TestUtils::ToVector( + cuda_stream, device_geometries.get_offsets().polygon_offsets.ps_num_rings); + ctx.prefix_sum2 = TestUtils::ToVector( + cuda_stream, device_geometries.get_offsets().polygon_offsets.ps_num_points); + ctx.points = TestUtils::ToVector(cuda_stream, device_geometries.get_points()); + ctx.mbrs = TestUtils::ToVector(cuda_stream, device_geometries.get_mbrs()); + cuda_stream.synchronize(); + PolygonArrayView<POINT_T, INDEX_T> poly_array_view( + ArrayView<INDEX_T>(ctx.prefix_sum1.data(), ctx.prefix_sum1.size()), + ArrayView<INDEX_T>(ctx.prefix_sum2.data(), ctx.prefix_sum2.size()), + ArrayView<POINT_T>(ctx.points.data(), ctx.points.size()), + ArrayView<box_t>(ctx.mbrs.data(), ctx.mbrs.size())); + poly = poly_array_view[0]; +} + +template <typename POINT_T, typename INDEX_T> +void ParseWKTMultiPolygon(Context<POINT_T, INDEX_T>& ctx, const char* wkt, + MultiPolygon<POINT_T, INDEX_T>& poly) { + nanoarrow::UniqueArrayStream stream; + ArrayStreamFromWKT({{wkt}}, GEOARROW_TYPE_WKB, stream.get()); + nanoarrow::UniqueArray array; + ArrowError error; + ArrowErrorSet(&error, ""); + + ASSERT_EQ(ArrowArrayStreamGetNext(stream.get(), array.get(), &error), NANOARROW_OK); + loader_t loader; + auto cuda_stream = rmm::cuda_stream_default; + + loader.Init(); + loader.Parse(cuda_stream, array.get(), 0, array->length); + auto device_geometries = loader.Finish(cuda_stream); + ctx.prefix_sum1 = TestUtils::ToVector( + cuda_stream, device_geometries.get_offsets().multi_polygon_offsets.ps_num_parts); + ctx.prefix_sum2 = TestUtils::ToVector( + cuda_stream, device_geometries.get_offsets().multi_polygon_offsets.ps_num_rings); + ctx.prefix_sum3 = TestUtils::ToVector( + cuda_stream, device_geometries.get_offsets().multi_polygon_offsets.ps_num_points); + ctx.points = TestUtils::ToVector(cuda_stream, device_geometries.get_points()); + ctx.mbrs = TestUtils::ToVector(cuda_stream, device_geometries.get_mbrs()); + cuda_stream.synchronize(); + MultiPolygonArrayView<POINT_T, INDEX_T> poly_array_view( + ArrayView<INDEX_T>(ctx.prefix_sum1.data(), ctx.prefix_sum1.size()), + ArrayView<INDEX_T>(ctx.prefix_sum2.data(), ctx.prefix_sum2.size()), + ArrayView<INDEX_T>(ctx.prefix_sum3.data(), ctx.prefix_sum3.size()), + ArrayView<POINT_T>(ctx.points.data(), ctx.points.size()), + ArrayView<box_t>(ctx.mbrs.data(), ctx.mbrs.size())); + poly = poly_array_view[0]; +} + +template <typename GEOMETRY1_T, typename GEOMETRY2_T> +void TestRelate(const char* wkt1, const char* wkt2, const GEOMETRY1_T& g1, + const GEOMETRY2_T& g2) { + WKTReader r; + auto a = r.read(wkt1); + auto b = r.read(wkt2); + + RelateMatrixPredicate pred; + RelateNG::relate(a.get(), b.get(), pred); + std::string actualVal = pred.getIM()->toString(); + + int val = relate(g1, g2); + char res[10]; + IM__ToString(val, res); + ASSERT_STREQ(actualVal.c_str(), res); +} + +TEST(RelateTest, PointPointDisjoint) { + point_t p1, p2; + + std::string wkt1 = "POINT (0 0)"; + std::string wkt2 = "POINT (1 1)"; + ParseWKTPoint(wkt1.c_str(), p1); + ParseWKTPoint(wkt2.c_str(), p2); + TestRelate(wkt1.c_str(), wkt2.c_str(), p1, p2); +} + +TEST(RelateTest, MultiPointMultiPointContained) { + MultiPoint<point_t> m1, m2; + std::string wkt1 = "MULTIPOINT (0 0, 1 1, 2 2)"; + std::string wkt2 = "MULTIPOINT (1 1, 2 2)"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTMultiPoint(ctx1, wkt1.c_str(), m1); + ParseWKTMultiPoint(ctx2, wkt2.c_str(), m2); + + TestRelate(wkt1.c_str(), wkt2.c_str(), m1, m2); +} + +TEST(RelateTest, MultiPointMultiPointEqual) { + using point_t = Point<double, 2>; + MultiPoint<point_t> m1, m2; + std::string wkt1 = "MULTIPOINT (0 0, 1 1, 2 2)"; + std::string wkt2 = "MULTIPOINT (0 0, 1 1, 2 2)"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTMultiPoint(ctx1, wkt1.c_str(), m1); + ParseWKTMultiPoint(ctx2, wkt2.c_str(), m2); + + TestRelate(wkt1.c_str(), wkt2.c_str(), m1, m2); +} + +TEST(RelateTest, MultiPointMultiPointValidateRelatePP_13) { + MultiPoint<point_t> m1, m2; + std::string wkt1 = "MULTIPOINT ((80 70), (140 120), (20 20), (200 170))"; + std::string wkt2 = "MULTIPOINT ((80 70), (140 120), (80 170), (200 80))"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTMultiPoint(ctx1, wkt1.c_str(), m1); + ParseWKTMultiPoint(ctx2, wkt2.c_str(), m2); + + TestRelate(wkt1.c_str(), wkt2.c_str(), m1, m2); +} + +TEST(RelateTest, LineStringMultiPointContains) { + LineString<point_t> ls1; + MultiPoint<point_t> m2; + std::string wkt1 = "LINESTRING (0 0, 1 1, 2 2)"; + std::string wkt2 = "MULTIPOINT (0 0, 1 1, 2 2)"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTLineString(ctx1, wkt1.c_str(), ls1); + ParseWKTMultiPoint(ctx2, wkt2.c_str(), m2); + TestRelate(wkt1.c_str(), wkt2.c_str(), ls1, m2); +} + +TEST(RelateTest, LineStringMultiPointOverlaps) { + LineString<point_t> ls1; + MultiPoint<point_t> m2; + std::string wkt1 = "LINESTRING (0 0, 1 1)"; + std::string wkt2 = "MULTIPOINT (0 0, 1 1, 2 2)"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTLineString(ctx1, wkt1.c_str(), ls1); + ParseWKTMultiPoint(ctx2, wkt2.c_str(), m2); + TestRelate(wkt1.c_str(), wkt2.c_str(), ls1, m2); +} + +TEST(RelateTest, ZeroLengthLinePoint) { + LineString<point_t> ls1; + point_t p2; + std::string wkt1 = "LINESTRING (0 0, 0 0)"; + std::string wkt2 = "POINT (0 0)"; + Context<point_t, index_t> ctx1; + + ParseWKTLineString(ctx1, wkt1.c_str(), ls1); + ParseWKTPoint(wkt2.c_str(), p2); + TestRelate(wkt1.c_str(), wkt2.c_str(), ls1, p2); +} + +TEST(RelateTest, ZeroLengthLineLine) { + LineString<point_t> ls1; + LineString<point_t> ls2; + std::string wkt1 = "LINESTRING (10 10, 10 10, 10 10)"; + std::string wkt2 = "LINESTRING (10 10, 10 10)"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTLineString(ctx1, wkt1.c_str(), ls1); + ParseWKTLineString(ctx2, wkt2.c_str(), ls2); + TestRelate(wkt1.c_str(), wkt2.c_str(), ls1, ls2); +} + +TEST(RelateTest, NonZeroLengthLinePoint) { + LineString<point_t> ls1; + point_t p2; + std::string wkt1 = "LINESTRING (0 0, 0 0, 9 9)"; + std::string wkt2 = "POINT (1 1)"; + Context<point_t, index_t> ctx1; + + ParseWKTLineString(ctx1, wkt1.c_str(), ls1); + ParseWKTPoint(wkt2.c_str(), p2); + TestRelate(wkt1.c_str(), wkt2.c_str(), ls1, p2); +} + +TEST(RelateTest, LinePointIntAndExt) { + MultiPoint<point_t> m1; + LineString<point_t> ls2; + std::string wkt1 = "MULTIPOINT ((60 60), (100 100))"; + std::string wkt2 = "LINESTRING (40 40, 80 80)"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTMultiPoint(ctx1, wkt1.c_str(), m1); + ParseWKTLineString(ctx2, wkt2.c_str(), ls2); + TestRelate(wkt1.c_str(), wkt2.c_str(), m1, ls2); +} + +TEST(RelateTest, LinesCrossProper) { + LineString<point_t> ls1; + LineString<point_t> ls2; + std::string wkt1 = "LINESTRING (0 0, 9 9)"; + std::string wkt2 = "LINESTRING (0 9, 9 0)"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTLineString(ctx1, wkt1.c_str(), ls1); + ParseWKTLineString(ctx2, wkt2.c_str(), ls2); + TestRelate(wkt1.c_str(), wkt2.c_str(), ls1, ls2); +} + +TEST(RelateTest, LinesOverlap) { + LineString<point_t> ls1; + LineString<point_t> ls2; + std::string wkt1 = "LINESTRING (0 0, 5 5)"; + std::string wkt2 = "LINESTRING (3 3, 9 9)"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTLineString(ctx1, wkt1.c_str(), ls1); + ParseWKTLineString(ctx2, wkt2.c_str(), ls2); + TestRelate(wkt1.c_str(), wkt2.c_str(), ls1, ls2); +} + +TEST(RelateTest, LinesCrossVertex) { + LineString<point_t> ls1; + LineString<point_t> ls2; + std::string wkt1 = "LINESTRING (0 0, 8 8)"; + std::string wkt2 = "LINESTRING (0 8, 4 4, 8 0)"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTLineString(ctx1, wkt1.c_str(), ls1); + ParseWKTLineString(ctx2, wkt2.c_str(), ls2); + TestRelate(wkt1.c_str(), wkt2.c_str(), ls1, ls2); +} + +TEST(RelateTest, LinesTouchVertex) { + LineString<point_t> ls1; + LineString<point_t> ls2; + std::string wkt1 = "LINESTRING (0 0, 8 0)"; + std::string wkt2 = "LINESTRING (0 8, 4 0, 8 8)"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTLineString(ctx1, wkt1.c_str(), ls1); + ParseWKTLineString(ctx2, wkt2.c_str(), ls2); + TestRelate(wkt1.c_str(), wkt2.c_str(), ls1, ls2); +} + +TEST(RelateTest, LinesDisjointByEnvelope) { + LineString<point_t> ls1; + LineString<point_t> ls2; + std::string wkt1 = "LINESTRING (0 0, 9 9)"; + std::string wkt2 = "LINESTRING (10 19, 19 10)"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTLineString(ctx1, wkt1.c_str(), ls1); + ParseWKTLineString(ctx2, wkt2.c_str(), ls2); + TestRelate(wkt1.c_str(), wkt2.c_str(), ls1, ls2); +} + +TEST(RelateTest, LinesDisjoint) { + LineString<point_t> ls1; + LineString<point_t> ls2; + std::string wkt1 = "LINESTRING (0 0, 9 9)"; + std::string wkt2 = "LINESTRING (4 2, 8 6)"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTLineString(ctx1, wkt1.c_str(), ls1); + ParseWKTLineString(ctx2, wkt2.c_str(), ls2); + TestRelate(wkt1.c_str(), wkt2.c_str(), ls1, ls2); +} + +// FIXME: wkt1 is a closed polyline, which has no boundary according to JTS's +// Mod2BoundaryNodeRule We have to implement a similar rule in gpuspatial to handle this +// case correctly TEST(RelateTest, LinesClosedEmpty) { +// MultiLineString<point_t, index_t> m_ls1; +// LineString<point_t> ls2; +// std::string wkt1 = "MULTILINESTRING ((0 0, 0 1), (0 1, 1 1, 1 0, 0 0))"; +// std::string wkt2 = "LINESTRING EMPTY"; +// Context<point_t, index_t> ctx1, ctx2; +// +// ParseWKTMultiLineString(ctx1, wkt1.c_str(), m_ls1); +// ParseWKTLineString(ctx2, wkt2.c_str(), ls2); +// TestRelate(wkt1.c_str(), wkt2.c_str(), m_ls1, ls2); +// } + +TEST(RelateTest, LinesRingTouchAtNode) { + LineString<point_t> ls1; + LineString<point_t> ls2; + std::string wkt1 = "LINESTRING (5 5, 1 8, 1 1, 5 5)"; + std::string wkt2 = "LINESTRING (5 5, 9 5)"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTLineString(ctx1, wkt1.c_str(), ls1); + ParseWKTLineString(ctx2, wkt2.c_str(), ls2); + TestRelate(wkt1.c_str(), wkt2.c_str(), ls1, ls2); +} + +TEST(RelateTest, LinesTouchAtBdy) { + LineString<point_t> ls1; + LineString<point_t> ls2; + std::string wkt1 = "LINESTRING (5 5, 1 8)"; + std::string wkt2 = "LINESTRING (5 5, 9 5)"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTLineString(ctx1, wkt1.c_str(), ls1); + ParseWKTLineString(ctx2, wkt2.c_str(), ls2); + TestRelate(wkt1.c_str(), wkt2.c_str(), ls1, ls2); +} + +TEST(RelateTest, LinesOverlapWithDisjointLine) { + LineString<point_t> ls1; + MultiLineString<point_t, index_t> m_ls2; + std::string wkt1 = "LINESTRING (1 1, 9 9)"; + std::string wkt2 = "MULTILINESTRING ((2 2, 8 8), (6 2, 8 4))"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTLineString(ctx1, wkt1.c_str(), ls1); + ParseWKTMultiLineString(ctx2, wkt2.c_str(), m_ls2); + TestRelate(wkt1.c_str(), wkt2.c_str(), ls1, m_ls2); +} + +TEST(RelateTest, LinesDisjointOverlappingEnvelopes) { + LineString<point_t> ls1; + LineString<point_t> ls2; + std::string wkt1 = "LINESTRING (60 0, 20 80, 100 80, 80 120, 40 140)"; + std::string wkt2 = "LINESTRING (60 40, 140 40, 140 160, 0 160)"; + Context<point_t, index_t> ctx1, ctx2; + + ParseWKTLineString(ctx1, wkt1.c_str(), ls1); + ParseWKTLineString(ctx2, wkt2.c_str(), ls2); + TestRelate(wkt1.c_str(), wkt2.c_str(), ls1, ls2); +} + +/** + * FIXME: + * Case from https://github.com/locationtech/jts/issues/270 + * Strictly, the lines cross, since their interiors intersect + * according to the Orientation predicate. + * However, the computation of the intersection point is + * non-robust, and reports it as being equal to the endpoint + * POINT (-10 0.0000000000000012) + * For consistency the relate algorithm uses the intersection node topology. + */ +// TEST(RelateTest, LineStringLineString10) { +// LineString<point_t> ls1; +// LineString<point_t> ls2; +// std::string wkt1 = "LINESTRING (0 0, -10 0.0000000000000012)"; +// std::string wkt2 = "LINESTRING (-9.999143275740073 -0.1308959557133398, -10 +// 0.0000000000001054)"; Context<point_t, index_t> ctx1, ctx2; +// +// ParseWKTLineString(ctx1, wkt1.c_str(), ls1); +// ParseWKTLineString(ctx2, wkt2.c_str(), ls2); +// TestRelate(wkt1.c_str(), wkt2.c_str(), ls1, ls2); +// } Review Comment: https://github.com/apache/sedona-db/issues/415 -- This is an automated message from the Apache Git Service. To respond to the message, please log on to GitHub and use the URL above to go to the specific comment. To unsubscribe, e-mail: [email protected] For queries about this service, please contact Infrastructure at: [email protected]
