OIiveirra commented on code in PR #69032:
URL: https://github.com/apache/doris/pull/69032#discussion_r4374955557
##########
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:
Fixed in 5dcdbae409e. Shared-arc noding now uses exact original-degree
coplanarity and longitude ordering, preserving equivalent subdivisions without
merging nearby curves. Added rotated, date-line, near-antipodal and tiny-offset
coverage; 110 geo/WKB UTs and both targeted regression suites pass.
##########
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:
Fixed in 5dcdbae409e. Meridian coverage is merged before noding, and atomic
edges are deduplicated while preserving original line boundaries. Added
query-accounted scratch, increasing-size retrace/dense-line storage-budget
tests, and SQL coverage; 110 geo/WKB UTs and both targeted regression suites
pass.
--
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]
