github-actions[bot] commented on code in PR #69032:
URL: https://github.com/apache/doris/pull/69032#discussion_r4379422808
##########
be/src/exprs/function/geo/functions_geo.cpp:
##########
@@ -18,32 +18,1175 @@
#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/s2contains_point_query.h>
+#include <s2/s2earth.h>
+#include <s2/s2lax_polygon_shape.h>
+#include <s2/s2loop.h>
+#include <s2/s2loop_measures.h>
#include <s2/s2point.h>
+#include <s2/s2point_index.h>
+#include <s2/s2point_vector_shape.h>
+#include <s2/s2pointutil.h>
+#include <s2/s2polygon.h>
+#include <s2/s2polyline.h>
+#include <s2/s2predicates.h>
+#include <s2/s2shape_measures.h>
#include <algorithm>
+#include <array>
+#include <bit>
#include <boost/iterator/iterator_facade.hpp>
+#include <boost/multiprecision/cpp_int.hpp>
+#include <cmath>
+#include <cstdint>
+#include <limits>
+#include <map>
+#include <numbers>
+#include <optional>
+#include <set>
#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/custom_allocator.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;
+}
+
+// Error-free TwoSum preserves tiny differences from a half-turn (e.g. -1e-15
+ 180).
+static std::pair<double, double> geography_degree_sum(double a, double b) {
+ const double high = a + b;
+ const double offset = high - a;
+ return {high, (a - (high - offset)) + (b - offset)};
+}
+
+// Reduce the exact degree sum before converting to radians. In particular,
+// 180 + (-1e-300) must retain its small sine instead of becoming sin(pi).
+static std::pair<long double, long double> geography_degree_sincos(double a,
double b = 0,
+ bool half =
false) {
+ const auto [high, low] = geography_degree_sum(a, b);
+ int quadrant;
+ const long double reduced =
+ std::remquo(static_cast<long double>(high), half ? 180.0L : 90.0L,
&quadrant) +
+ static_cast<long double>(low);
+ const auto radians = reduced * (std::numbers::pi_v<long double> / 180) /
(half ? 2 : 1);
+ const auto sine = std::sin(radians);
+ const auto cosine = std::cos(radians);
+ switch ((quadrant % 4 + 4) % 4) {
+ case 0:
+ return {sine, cosine};
+ case 1:
+ return {cosine, -sine};
+ case 2:
+ return {-sine, -cosine};
+ default:
+ return {-cosine, sine};
+ }
+}
+
+// Unlike cross(ToPoint(a), ToPoint(b)), these product identities do not
subtract
+// rounded nearly opposite vectors. Longitude-midpoint coordinates give
+// (-sin(dLon/2)*sin(latA+latB), cos(dLon/2)*sin(latA-latB),
+// cos(latA)*cos(latB)*sin(dLon)) for the oriented normal.
+struct GeographyRawArc {
+ S2Point normal;
+ S2Point midpoint;
+
+ GeographyRawArc(const GeoCoordinate& a, const GeoCoordinate& b) {
+ // Wide exponents retain products of the smallest binary64 degree
residuals.
+ static_assert(std::numeric_limits<long double>::digits >= 64 &&
+ std::numeric_limits<long double>::min_exponent <= -4096);
+ const auto [sd, cd] = geography_degree_sincos(b.x, -a.x, true);
+ const auto [sm, cm] = geography_degree_sincos(a.x, b.x, true);
+ const auto rotate_unit = [&](long double x, long double y, long double
z) {
+ const auto norm = std::sqrt(x * x + y * y + z * z);
+ DCHECK_GT(norm, 0);
+ return S2Point(static_cast<double>((cm * x - sm * y) / norm),
+ static_cast<double>((sm * x + cm * y) / norm),
+ static_cast<double>(z / norm));
+ };
+ normal = rotate_unit(-sd * geography_degree_sincos(a.y, b.y).first,
+ cd * geography_degree_sincos(a.y, -b.y).first,
+ geography_degree_sincos(a.y).second *
+ geography_degree_sincos(b.y).second *
+ geography_degree_sincos(b.x, -a.x).first);
+ const auto [sl, cl] = geography_degree_sincos(a.y, b.y, true);
+ const auto [st, ct] = geography_degree_sincos(b.y, -a.y, true);
+ midpoint = rotate_unit(cd * cl * ct, -sd * sl * st, sl * ct);
+ }
+};
+
+#ifdef BE_TEST
+static thread_local size_t geography_exact_proofs = 0;
+static thread_local size_t geography_candidate_visits = 0;
+#endif
+
+// Test coplanarity of the original degree ordinates, not independently
rounded S2 vectors.
+// For example (45,-45), (90,0), (135,45) share an arc although ToPoint's
determinant is nonzero.
+static bool geography_common_great_circle(const GeoCoordinate& a, const
GeoCoordinate& b,
+ const GeoCoordinate& c) {
+#ifdef BE_TEST
+ ++geography_exact_proofs;
+#endif
+ using boost::multiprecision::cpp_int;
+ const std::array<double, 6> ordinates {a.y, a.x, b.y, b.x, c.y, c.x};
+ std::array<uint64_t, 6> significands;
+ std::array<int, 6> exponents;
+ int denominator_bits = 0;
+ for (size_t i = 0; i < ordinates.size(); ++i) {
+ const auto bits = std::bit_cast<uint64_t>(ordinates[i]);
+ const int biased = (bits >> 52) & 0x7ff;
+ uint64_t significand = bits & ((uint64_t {1} << 52) - 1);
+ int exponent = biased ? biased - 1075 : -1074;
+ if (biased) {
+ significand |= uint64_t {1} << 52;
+ }
+ if (significand != 0) {
+ const int zeros = std::countr_zero(significand);
+ significand >>= zeros;
+ exponent += zeros;
+ denominator_bits = std::max(denominator_bits, -exponent);
+ } else {
+ exponent = 0;
+ }
+ significands[i] = significand;
+ exponents[i] = exponent;
+ }
+ std::array<cpp_int, 6> angles;
+ for (size_t i = 0; i < angles.size(); ++i) {
+ angles[i] = cpp_int(significands[i]) << (denominator_bits +
exponents[i]);
+ if (std::signbit(ordinates[i])) {
+ angles[i] = -angles[i];
+ }
+ }
+ // Binary-rational degrees have denominator 2^q. Exp(i*angle) is a power
of a
+ // root of unity of order 360*2^q. Its minimal polynomial is Phi_90(t),
with
+ // t = z^(2^(q+2)), and Phi_90(t) = t^24+t^21-t^15-t^12-t^9+t^3+1.
+ // Expand cos(lat_i)cos(lat_j)sin(lat_k)sin(lon_j-lon_i), cyclically, into
+ // 48 integer-coefficient powers. Equal residues modulo 2^(q+2) share one
+ // polynomial; reduction by Phi_90 proves zero exactly, including tiny
nonzero controls.
+ const int shift = denominator_bits + 2;
+ const cpp_int basis = cpp_int(1) << shift;
+ const cpp_int period = 90 * basis;
+ std::map<cpp_int, std::array<int, 90>> polynomials;
+ for (int i = 0; i < 3; ++i) {
+ const int j = (i + 1) % 3;
+ const int k = (i + 2) % 3;
+ for (int mask = 0; mask < 16; ++mask) {
+ std::array<int, 4> signs;
+ for (int bit = 0; bit < 4; ++bit) {
+ signs[bit] = (mask & (1 << bit)) ? 1 : -1;
+ }
+ cpp_int power = signs[0] * angles[2 * i] + signs[1] * angles[2 *
j] +
+ signs[2] * angles[2 * k] +
+ signs[3] * (angles[2 * j + 1] - angles[2 * i + 1]);
+ power %= period;
+ if (power < 0) {
+ power += period;
+ }
+ const int quotient = (power >> shift).convert_to<int>();
+ const cpp_int residue = power - (cpp_int(quotient) << shift);
+ polynomials[residue][quotient] += signs[2] * signs[3];
+ }
+ }
+ for (auto& [residue, coefficients] : polynomials) {
+ for (int i = 89; i >= 24; --i) {
+ const int coefficient = coefficients[i];
+ coefficients[i] = 0;
+ coefficients[i - 3] -= coefficient;
+ coefficients[i - 9] += coefficient;
+ coefficients[i - 12] += coefficient;
+ coefficients[i - 15] += coefficient;
+ coefficients[i - 21] -= coefficient;
+ coefficients[i - 24] -= coefficient;
+ }
+ if (std::ranges::any_of(coefficients, [](int value) { return value !=
0; })) {
+ return false;
+ }
+ }
+ return true;
+}
+
+// Keep each atomic segment in its own shape. Pinned S2 ProcessEdge1 uses a
chain's
+// global start as a local offset in one endpoint path, so multiple chains per
line
+// shape are unsafe. Closed AA shapes restore internal vertices in the OPEN
model.
+class GeographyLineShape final : public S2Shape {
+public:
+ GeographyLineShape(const S2Point& a, const S2Point& b) {
_edge.emplace_back(a, b); }
+ int num_edges() const override { return 1; }
+ Edge edge(int i) const override {
+ DCHECK_EQ(i, 0);
+ return _edge.front();
+ }
+ int dimension() const override { return 1; }
+ ReferencePoint GetReferencePoint() const override { return
ReferencePoint::Contained(false); }
+ int num_chains() const override { return 1; }
+ Chain chain(int i) const override {
+ DCHECK_EQ(i, 0);
+ return {0, 1};
+ }
+ Edge chain_edge(int i, int j) const override {
+ DCHECK_EQ(i, 0);
+ return edge(j);
+ }
+ ChainPosition chain_position(int i) const override {
+ DCHECK_EQ(i, 0);
+ return {0, 0};
+ }
+
+private:
+ // Owning storage is checked/accounted by the current query's allocator.
+ DorisVector<Edge> _edge;
+};
+
+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);
+ DorisMap<double, DorisVector<S2Point>> longitude_vertices;
+ 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());
+ if (!_coordinates.emplace(point, coordinate).second) {
Review Comment:
[P2] Preserve raw endpoint identity when rounded points collide
For default GEOGRAPHY WKB, let P=`POINT(0 58)` and L=`LINESTRING(0
58.00000000000001,180 -58)`. The pinned [S1Angle degree
conversion](https://github.com/google/s2geometry/blob/v0.10.0/src/s2/s1angle.h#L290-L292)
maps both distinct latitudes to the same binary64 radians
`0x1.03258031ade7cp+0`, so their S2Points coincide. This constructor processes
P first and discards L's original start coordinate. `add_line -> meridian_edge`
then recovers endpoints `(0,58)` and `(-180,-58)`, triggering the
exact-antipodal exception in `ST_Distance(P,L)` and `ST_DWithin(P,L,1)`.
Reversing operands retains L's original latitude and avoids that guard. L
itself is non-antipodal and passes the [pinned polyline
validation](https://github.com/google/s2geometry/blob/v0.10.0/src/s2/s2polyline.cc#L148-L169).
This is source/arithmetic evidence, not a Doris runtime test. It is a
distinct lossy-identity failure upstream of existing distance-plane finding
4376498668. Preserve edge-local canonical raw coordinates independently of the
rounded S2Point key, and cover colliding latitudes in both operand orders while
retaining pole/date-line aliases.
##########
be/src/exprs/function/geo/functions_geo.cpp:
##########
@@ -18,32 +18,1175 @@
#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/s2contains_point_query.h>
+#include <s2/s2earth.h>
+#include <s2/s2lax_polygon_shape.h>
+#include <s2/s2loop.h>
+#include <s2/s2loop_measures.h>
#include <s2/s2point.h>
+#include <s2/s2point_index.h>
+#include <s2/s2point_vector_shape.h>
+#include <s2/s2pointutil.h>
+#include <s2/s2polygon.h>
+#include <s2/s2polyline.h>
+#include <s2/s2predicates.h>
+#include <s2/s2shape_measures.h>
#include <algorithm>
+#include <array>
+#include <bit>
#include <boost/iterator/iterator_facade.hpp>
+#include <boost/multiprecision/cpp_int.hpp>
+#include <cmath>
+#include <cstdint>
+#include <limits>
+#include <map>
+#include <numbers>
+#include <optional>
+#include <set>
#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/custom_allocator.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);
Review Comment:
[P2] Stabilize typed polygon holes before legacy admission
Import this default GEOGRAPHY WKB polygon:
```text
POLYGON((-44.99999999999997 -45,134.99999999999997 45,45
45,-44.99999999999997 -45),(44.9999999 0.0000009,45.0000001 0.0000009,45
0.0000011,44.9999999 0.0000009))
```
Its small hole lies strictly inside the raw spherical shell. This decode
still reaches `to_s2polygon` (`geo_types.cpp:384`), whose rounded-loop
hole-containment check runs before
`GeographyEdgeNodes`/`geography_shape_index`. Modeling the pinned [ToPoint
conversion](https://github.com/google/s2geometry/blob/v0.10.0/src/s2/s2latlng.cc#L32-L38)
with exact determinants of its binary64 outputs gives negative
first-shell-edge signs for all three hole vertices (about `-2.77556e-17`); the
other two sides are positive. Thus
[S2Loop::Contains](https://github.com/google/s2geometry/blob/v0.10.0/src/s2/s2loop.cc#L955-L988)
rejects the hole, and typed area/predicates report an unsupported shape.
Inserting the exact raw midpoint `(45,0)` on that shell edge preserves the
polygon but makes both split-edge signs positive (minimum `9.87307e-9`) and
admission succeeds.
This is source/arithmetic evidence, not a runtime test. Unlike existing
interior/area findings 4377327193/4378237155, this upstream consumer rejects
the whole polygon before their repaired calculations run; current hole tests at
latitudes 20-22 miss it. Stabilize typed rings before hole admission and cover
near-edge holes with equivalent subdivision/reversal.
--
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]
