github-actions[bot] commented on code in PR #69032:
URL: https://github.com/apache/doris/pull/69032#discussion_r4374330645
##########
be/src/exprs/function/geo/functions_geo.cpp:
##########
@@ -18,32 +18,634 @@
#include "exprs/function/geo/functions_geo.h"
#include <glog/logging.h>
+#include <s2/mutable_s2shape_index.h>
+#include <s2/s2boolean_operation.h>
+#include <s2/s2builderutil_s2polyline_vector_layer.h>
+#include <s2/s2closest_edge_query.h>
+#include <s2/s2closest_point_query.h>
+#include <s2/s2earth.h>
+#include <s2/s2lax_polygon_shape.h>
+#include <s2/s2lax_polyline_shape.h>
+#include <s2/s2loop.h>
#include <s2/s2point.h>
+#include <s2/s2point_index.h>
+#include <s2/s2point_vector_shape.h>
+#include <s2/s2polygon.h>
+#include <s2/s2polyline.h>
+#include <s2/s2predicates.h>
#include <algorithm>
+#include <array>
#include <boost/iterator/iterator_facade.hpp>
+#include <cmath>
+#include <cstdint>
+#include <limits>
+#include <map>
+#include <optional>
#include <utility>
+#include "common/cast_set.h"
#include "common/compiler_util.h"
#include "core/assert_cast.h"
#include "core/block/block.h"
#include "core/block/column_with_type_and_name.h"
#include "core/column/column.h"
#include "core/column/column_array.h"
#include "core/column/column_execute_util.h"
+#include "core/column/column_map.h"
#include "core/column/column_nullable.h"
+#include "core/column/column_struct.h"
+#include "core/column/column_varbinary.h"
#include "core/data_type/data_type_array.h"
+#include "core/data_type/data_type_map.h"
#include "core/data_type/data_type_nullable.h"
#include "core/data_type/data_type_number.h"
+#include "core/data_type/data_type_spatial.h"
#include "core/data_type/data_type_string.h"
+#include "core/data_type/data_type_struct.h"
+#include "core/data_type/data_type_varbinary.h"
#include "core/data_type/define_primitive_type.h"
#include "core/string_ref.h"
#include "exprs/function/geo/geo_common.h"
#include "exprs/function/geo/geo_types.h"
+#include "exprs/function/geo/wkb_parse.h"
#include "exprs/function/simple_function_factory.h"
+#include "exprs/function/string_hex_util.h"
namespace doris {
+static bool is_spatial_type(const DataTypePtr& type) {
+ const auto primitive_type = remove_nullable(type)->get_primitive_type();
+ return primitive_type == TYPE_GEOMETRY || primitive_type == TYPE_GEOGRAPHY;
+}
+
+static Status validate_geography_semantics(const DataTypePtr& type, const
char* function_name) {
+ const auto& nested_type = remove_nullable(type);
+ if (!is_spatial_type(nested_type)) {
+ return Status::OK();
+ }
+
+ const auto* spatial_type = dynamic_cast<const
DataTypeSpatial*>(nested_type.get());
+ DCHECK(spatial_type != nullptr);
+ if (spatial_type->get_primitive_type() == TYPE_GEOGRAPHY &&
+ spatial_type->crs() == "OGC:CRS84" && spatial_type->algorithm() ==
"spherical") {
+ return Status::OK();
+ }
+ return Status::NotSupported(
+ "Function {} requires GEOGRAPHY(OGC:CRS84, spherical) for spatial
inputs",
+ function_name);
+}
+
+static std::unique_ptr<GeoShape> decode_geo_shape(StringRef value, const
DataTypePtr& type) {
+ if (!is_spatial_type(type)) {
+ return GeoShape::from_encoded(value.data, value.size);
+ }
+
+ GeoParseStatus status;
+ auto shape = GeoShape::from_wkb_bytes(value.data, value.size, status);
+ return status == GeoParseStatus::GEO_PARSE_OK ? std::move(shape) : nullptr;
+}
+
+static bool any_geo_argument_is_null(const Block& block, const ColumnNumbers&
arguments,
+ size_t row) {
+ return std::ranges::any_of(arguments, [&](auto argument) {
+ return block.get_by_position(argument).column->is_null_at(row);
+ });
+}
+
+static Status decode_calculation_shape(StringRef value, const DataTypePtr&
type,
+ std::unique_ptr<GeoShape>* shape,
+ const char* function_name) {
+ *shape = decode_geo_shape(value, type);
+ // POLYGON EMPTY decodes successfully, but distance helpers retain
infinity/max sentinels.
+ // Reject it for typed calculations just like empty points/lines, leaving
raw WKB untouched.
+ if (is_spatial_type(type) &&
+ (*shape == nullptr ||
+ ((*shape)->type() == GeoShapeType::GEO_SHAPE_POLYGON &&
+ assert_cast<const GeoPolygon&>(**shape).polygon()->num_loops() ==
0))) {
+ return Status::NotSupported("Function {} does not support this
Geography WKB shape",
+ function_name);
+ }
+ return Status::OK();
+}
+
+// Calculations use canonical positions; raw spatial accessors still return
the original WKB.
+static S2Point canonical_geography_point(const S2Point& point) {
+ const auto coordinate = S2LatLng(point);
+ const double latitude = coordinate.lat().radians();
+ // All longitudes at a pole and +/-180 at the date line denote the same
position.
+ if (std::abs(latitude) == M_PI_2) {
+ return {0, 0, latitude > 0 ? 1.0 : -1.0};
+ }
+ if (coordinate.lng().radians() == M_PI) {
+ // Flip only sin(longitude): rebuilding latitude can change its last
bit.
+ return {point.x(), -point.y(), point.z()};
+ }
+ return point;
+}
+
+// S2 predicates operate on exact rounded vectors. Splitting a meridian at a
new latitude
+// can otherwise turn one common edge into two almost-coplanar curves. Node
only edges whose
+// shared great circle is proved by original WKB ordinates; never snap nearby
parallel curves.
+class GeographyEdgeNodes {
+public:
+ GeographyEdgeNodes(StringRef left, StringRef right) {
+ GeoCoordinateList coordinates;
+ const auto left_status = WkbParse::xy_coordinates(left.data,
left.size, &coordinates);
+ const auto right_status = WkbParse::xy_coordinates(right.data,
right.size, &coordinates);
+ DCHECK_EQ(left_status, GeoParseStatus::GEO_PARSE_OK);
+ DCHECK_EQ(right_status, GeoParseStatus::GEO_PARSE_OK);
+ for (auto coordinate : coordinates.list) {
+ if (coordinate.x == 180) {
+ coordinate.x = -180;
+ }
+ const auto point = canonical_geography_point(
+ S2LatLng::FromDegrees(coordinate.y,
coordinate.x).ToPoint());
+ _coordinates.emplace(point, coordinate);
+ _vertex_index.Add(point, 0);
+ if (std::abs(coordinate.y) != 90) {
+ _meridians[meridian(coordinate.x)][coordinate.x < 0].push_back(
+ {coordinate.y, point});
+ }
+ }
+ for (auto& entry : _meridians) {
+ for (auto& side : entry.second) {
+ std::ranges::sort(side, {}, &Vertex::latitude);
+ const auto duplicate = std::ranges::unique(side, {},
&Vertex::point);
+ side.erase(duplicate.begin(), duplicate.end());
+ }
+ }
+ }
+
+ void append_edge(const S2Point& raw_start, const S2Point& raw_end,
+ std::vector<S2Point>* output) const {
+ const auto start = canonical_geography_point(raw_start);
+ const auto end = canonical_geography_point(raw_end);
+ output->push_back(start);
+ // Different raw longitudes at one pole denote the same calculation
point.
+ if (start == end) {
+ return;
+ }
+ const auto& a = _coordinates.at(start);
+ const auto& b = _coordinates.at(end);
+ const bool a_pole = std::abs(a.y) == 90;
+ const bool b_pole = std::abs(b.y) == 90;
+ const auto key = meridian(a_pole ? b.x : a.x);
+ const bool common_meridian = a_pole || b_pole || key == meridian(b.x);
+ // Exact antipodes have no unique edge; rounded S2 vectors can conceal
this ambiguity.
+ if (a.y == -b.y && common_meridian && ((a_pole && b_pole) || (a.x < 0)
!= (b.x < 0))) {
+ throw Exception(ErrorCode::NOT_IMPLEMENTED_ERROR,
+ "Geography calculations do not support antipodal
edges");
+ }
+ if (common_meridian) {
+ const bool a_side = (a_pole ? b.x : a.x) < 0;
+ const bool b_side = (b_pole ? a.x : b.x) < 0;
+ const auto& group = _meridians.at(key);
+ if (a_side == b_side) {
+ append_range(group[a_side], a.y, b.y, output);
+ } else {
+ // (10,80) -> (-170,80) passes through the north pole. Keep
latitude
+ // as-is on each side: computing 180-latitude would erase tiny
ordinates.
+ const double latitude = a.y > -b.y ? 90 : -90;
+ append_range(group[a_side], a.y, latitude, output);
+ output->emplace_back(0, 0, latitude > 0 ? 1 : -1);
+ append_range(group[b_side], latitude, b.y, output);
+ }
+ return;
+ }
+ // For other great circles, including the equator, node only exact
coplanar
+ // points on the edge. The query radius selects candidates, never
equality.
+ S2ClosestPointQuery<int> query(&_vertex_index);
+ query.mutable_options()->set_conservative_max_distance(
+ S1Angle::Radians(8 * std::numeric_limits<double>::epsilon()));
+ S2ClosestPointQuery<int>::EdgeTarget target(start, end);
+ std::vector<S2Point> on_edge;
+ for (const auto& candidate : query.FindClosestPoints(&target)) {
+ const auto& point = candidate.point();
+ if (point != start && point != end &&
+ s2pred::ExpensiveSign(start, end, point, false) == 0 &&
Review Comment:
[P2] Preserve non-meridian geodesic subdivisions in topology predicates
For default GEOGRAPHY values A=`LINESTRING(45 -45,135 45)` and
B=`LINESTRING(45 -45,90 0,135 45)` imported as WKB, `(90,0)` is exactly the
midpoint of the same 120-degree great-circle arc. However, `ToPoint()` rounds
the three vectors independently: their determinant is about `-3.52e-17`, so
this exact-coplanarity condition rejects the midpoint and A is never noded to
match B. `ST_Equals(A,B)`, mutual `ST_Contains`, and `ST_Within` consequently
return false. These are valid nonempty XY lines with no antipodal edge. S2's
[Boolean-output
path](https://github.com/google/s2geometry/blob/v0.10.0/src/s2/s2boolean_operation.cc#L2570-L2577)
computes the exact predicate before constructing its builder, so builder
snapping cannot repair this. The conversion-distance fallback below only
repairs Intersects/Disjoint; the prior point-membership and intersection
threads do not cover these new Equals/nonpoint-containment paths. Preserve
shared non-meridian arcs under coordinate conversion without
merging nearby distinct curves, and add rotated subdivision tests for
Equals/Contains/Within plus nearby controls. This result follows from static
code and exact arithmetic, not a Doris runtime execution.
##########
be/src/exprs/function/geo/functions_geo.cpp:
##########
@@ -18,32 +18,634 @@
#include "exprs/function/geo/functions_geo.h"
#include <glog/logging.h>
+#include <s2/mutable_s2shape_index.h>
+#include <s2/s2boolean_operation.h>
+#include <s2/s2builderutil_s2polyline_vector_layer.h>
+#include <s2/s2closest_edge_query.h>
+#include <s2/s2closest_point_query.h>
+#include <s2/s2earth.h>
+#include <s2/s2lax_polygon_shape.h>
+#include <s2/s2lax_polyline_shape.h>
+#include <s2/s2loop.h>
#include <s2/s2point.h>
+#include <s2/s2point_index.h>
+#include <s2/s2point_vector_shape.h>
+#include <s2/s2polygon.h>
+#include <s2/s2polyline.h>
+#include <s2/s2predicates.h>
#include <algorithm>
+#include <array>
#include <boost/iterator/iterator_facade.hpp>
+#include <cmath>
+#include <cstdint>
+#include <limits>
+#include <map>
+#include <optional>
#include <utility>
+#include "common/cast_set.h"
#include "common/compiler_util.h"
#include "core/assert_cast.h"
#include "core/block/block.h"
#include "core/block/column_with_type_and_name.h"
#include "core/column/column.h"
#include "core/column/column_array.h"
#include "core/column/column_execute_util.h"
+#include "core/column/column_map.h"
#include "core/column/column_nullable.h"
+#include "core/column/column_struct.h"
+#include "core/column/column_varbinary.h"
#include "core/data_type/data_type_array.h"
+#include "core/data_type/data_type_map.h"
#include "core/data_type/data_type_nullable.h"
#include "core/data_type/data_type_number.h"
+#include "core/data_type/data_type_spatial.h"
#include "core/data_type/data_type_string.h"
+#include "core/data_type/data_type_struct.h"
+#include "core/data_type/data_type_varbinary.h"
#include "core/data_type/define_primitive_type.h"
#include "core/string_ref.h"
#include "exprs/function/geo/geo_common.h"
#include "exprs/function/geo/geo_types.h"
+#include "exprs/function/geo/wkb_parse.h"
#include "exprs/function/simple_function_factory.h"
+#include "exprs/function/string_hex_util.h"
namespace doris {
+static bool is_spatial_type(const DataTypePtr& type) {
+ const auto primitive_type = remove_nullable(type)->get_primitive_type();
+ return primitive_type == TYPE_GEOMETRY || primitive_type == TYPE_GEOGRAPHY;
+}
+
+static Status validate_geography_semantics(const DataTypePtr& type, const
char* function_name) {
+ const auto& nested_type = remove_nullable(type);
+ if (!is_spatial_type(nested_type)) {
+ return Status::OK();
+ }
+
+ const auto* spatial_type = dynamic_cast<const
DataTypeSpatial*>(nested_type.get());
+ DCHECK(spatial_type != nullptr);
+ if (spatial_type->get_primitive_type() == TYPE_GEOGRAPHY &&
+ spatial_type->crs() == "OGC:CRS84" && spatial_type->algorithm() ==
"spherical") {
+ return Status::OK();
+ }
+ return Status::NotSupported(
+ "Function {} requires GEOGRAPHY(OGC:CRS84, spherical) for spatial
inputs",
+ function_name);
+}
+
+static std::unique_ptr<GeoShape> decode_geo_shape(StringRef value, const
DataTypePtr& type) {
+ if (!is_spatial_type(type)) {
+ return GeoShape::from_encoded(value.data, value.size);
+ }
+
+ GeoParseStatus status;
+ auto shape = GeoShape::from_wkb_bytes(value.data, value.size, status);
+ return status == GeoParseStatus::GEO_PARSE_OK ? std::move(shape) : nullptr;
+}
+
+static bool any_geo_argument_is_null(const Block& block, const ColumnNumbers&
arguments,
+ size_t row) {
+ return std::ranges::any_of(arguments, [&](auto argument) {
+ return block.get_by_position(argument).column->is_null_at(row);
+ });
+}
+
+static Status decode_calculation_shape(StringRef value, const DataTypePtr&
type,
+ std::unique_ptr<GeoShape>* shape,
+ const char* function_name) {
+ *shape = decode_geo_shape(value, type);
+ // POLYGON EMPTY decodes successfully, but distance helpers retain
infinity/max sentinels.
+ // Reject it for typed calculations just like empty points/lines, leaving
raw WKB untouched.
+ if (is_spatial_type(type) &&
+ (*shape == nullptr ||
+ ((*shape)->type() == GeoShapeType::GEO_SHAPE_POLYGON &&
+ assert_cast<const GeoPolygon&>(**shape).polygon()->num_loops() ==
0))) {
+ return Status::NotSupported("Function {} does not support this
Geography WKB shape",
+ function_name);
+ }
+ return Status::OK();
+}
+
+// Calculations use canonical positions; raw spatial accessors still return
the original WKB.
+static S2Point canonical_geography_point(const S2Point& point) {
+ const auto coordinate = S2LatLng(point);
+ const double latitude = coordinate.lat().radians();
+ // All longitudes at a pole and +/-180 at the date line denote the same
position.
+ if (std::abs(latitude) == M_PI_2) {
+ return {0, 0, latitude > 0 ? 1.0 : -1.0};
+ }
+ if (coordinate.lng().radians() == M_PI) {
+ // Flip only sin(longitude): rebuilding latitude can change its last
bit.
+ return {point.x(), -point.y(), point.z()};
+ }
+ return point;
+}
+
+// S2 predicates operate on exact rounded vectors. Splitting a meridian at a
new latitude
+// can otherwise turn one common edge into two almost-coplanar curves. Node
only edges whose
+// shared great circle is proved by original WKB ordinates; never snap nearby
parallel curves.
+class GeographyEdgeNodes {
+public:
+ GeographyEdgeNodes(StringRef left, StringRef right) {
+ GeoCoordinateList coordinates;
+ const auto left_status = WkbParse::xy_coordinates(left.data,
left.size, &coordinates);
+ const auto right_status = WkbParse::xy_coordinates(right.data,
right.size, &coordinates);
+ DCHECK_EQ(left_status, GeoParseStatus::GEO_PARSE_OK);
+ DCHECK_EQ(right_status, GeoParseStatus::GEO_PARSE_OK);
+ for (auto coordinate : coordinates.list) {
+ if (coordinate.x == 180) {
+ coordinate.x = -180;
+ }
+ const auto point = canonical_geography_point(
+ S2LatLng::FromDegrees(coordinate.y,
coordinate.x).ToPoint());
+ _coordinates.emplace(point, coordinate);
+ _vertex_index.Add(point, 0);
+ if (std::abs(coordinate.y) != 90) {
+ _meridians[meridian(coordinate.x)][coordinate.x < 0].push_back(
+ {coordinate.y, point});
+ }
+ }
+ for (auto& entry : _meridians) {
+ for (auto& side : entry.second) {
+ std::ranges::sort(side, {}, &Vertex::latitude);
+ const auto duplicate = std::ranges::unique(side, {},
&Vertex::point);
+ side.erase(duplicate.begin(), duplicate.end());
+ }
+ }
+ }
+
+ void append_edge(const S2Point& raw_start, const S2Point& raw_end,
+ std::vector<S2Point>* output) const {
+ const auto start = canonical_geography_point(raw_start);
+ const auto end = canonical_geography_point(raw_end);
+ output->push_back(start);
+ // Different raw longitudes at one pole denote the same calculation
point.
+ if (start == end) {
+ return;
+ }
+ const auto& a = _coordinates.at(start);
+ const auto& b = _coordinates.at(end);
+ const bool a_pole = std::abs(a.y) == 90;
+ const bool b_pole = std::abs(b.y) == 90;
+ const auto key = meridian(a_pole ? b.x : a.x);
+ const bool common_meridian = a_pole || b_pole || key == meridian(b.x);
+ // Exact antipodes have no unique edge; rounded S2 vectors can conceal
this ambiguity.
+ if (a.y == -b.y && common_meridian && ((a_pole && b_pole) || (a.x < 0)
!= (b.x < 0))) {
+ throw Exception(ErrorCode::NOT_IMPLEMENTED_ERROR,
+ "Geography calculations do not support antipodal
edges");
+ }
+ if (common_meridian) {
+ const bool a_side = (a_pole ? b.x : a.x) < 0;
+ const bool b_side = (b_pole ? a.x : b.x) < 0;
+ const auto& group = _meridians.at(key);
+ if (a_side == b_side) {
+ append_range(group[a_side], a.y, b.y, output);
+ } else {
+ // (10,80) -> (-170,80) passes through the north pole. Keep
latitude
+ // as-is on each side: computing 180-latitude would erase tiny
ordinates.
+ const double latitude = a.y > -b.y ? 90 : -90;
+ append_range(group[a_side], a.y, latitude, output);
+ output->emplace_back(0, 0, latitude > 0 ? 1 : -1);
+ append_range(group[b_side], latitude, b.y, output);
+ }
+ return;
+ }
+ // For other great circles, including the equator, node only exact
coplanar
+ // points on the edge. The query radius selects candidates, never
equality.
+ S2ClosestPointQuery<int> query(&_vertex_index);
+ query.mutable_options()->set_conservative_max_distance(
+ S1Angle::Radians(8 * std::numeric_limits<double>::epsilon()));
+ S2ClosestPointQuery<int>::EdgeTarget target(start, end);
+ std::vector<S2Point> on_edge;
+ for (const auto& candidate : query.FindClosestPoints(&target)) {
+ const auto& point = candidate.point();
+ if (point != start && point != end &&
+ s2pred::ExpensiveSign(start, end, point, false) == 0 &&
+ s2pred::CompareEdgeDistance(point, start, end,
S1ChordAngle::Zero()) <= 0) {
+ on_edge.push_back(point);
+ }
+ }
+ // Rounded angles can tie near the middle of a long edge, reversing
adjacent nodes.
+ std::ranges::sort(on_edge, [&](const S2Point& a, const S2Point& b) {
+ return s2pred::CompareDistances(start, a, b) < 0;
+ });
+ for (const auto& point : on_edge) {
+ if (output->back() != point) {
+ output->push_back(point);
+ }
+ }
+ }
+
+private:
+ struct Vertex {
+ double latitude;
+ S2Point point;
+ };
+ static std::pair<double, double> meridian(double longitude) {
+ if (longitude >= 0) {
+ return {longitude, 0};
+ }
+ // Error-free TwoSum: -1e-15+180 and -2e-15+180 both round to 180.
+ // Retain the residual so distinct original meridians never become one.
+ const double high = longitude + 180;
+ const double offset = high - 180;
+ const double low = (180 - (high - offset)) + (longitude - offset);
+ return {high, low};
+ }
+ static void append_range(const std::vector<Vertex>& vertices, double from,
double to,
+ std::vector<S2Point>* output) {
+ if (from == to) {
+ return;
+ }
+ auto first = std::ranges::upper_bound(vertices, std::min(from, to),
{}, &Vertex::latitude);
+ auto last = std::ranges::lower_bound(vertices, std::max(from, to), {},
&Vertex::latitude);
+ if (from < to) {
+ for (auto it = first; it != last; ++it) {
+ output->push_back(it->point);
Review Comment:
[P2] Avoid quadratic scratch expansion when noding meridian edges
Every common-meridian edge copies all vertices inside its latitude range.
Let A have N edges alternating `(10,0)` and `(10,1)`, and B have M distinct
vertices `(10,i/(M+1))`, `i=1..M`. These are admitted default-Geography lines:
only adjacent duplicates are removed, and [S2
validation](https://github.com/google/s2geometry/blob/v0.10.0/src/s2/s2polyline.cc#L160-L182)
permits retracing. Each A edge appends all M B vertices, so
`geography_shape_index` materializes `N*(M+1)+1` points before even
`ST_Intersects(A,B)` can return true. At N=M=10,000, the two WKB values total
320,034 bytes but the first vector alone contains 100,010,001 S2Points (at
least 2.4 GB of payload); larger valid lines quickly exhaust memory.
Deduplicating the candidate list does not reduce these repeated per-edge
copies. Avoid the full cross-product materialization and bound/account scratch
before allocating it; add retraced-line versus dense-collinear-line complexity
coverage. Counts are derived statically; no O
OM experiment or benchmark was run.
--
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]
---------------------------------------------------------------------
To unsubscribe, e-mail: [email protected]
For additional commands, e-mail: [email protected]
