github-actions[bot] commented on code in PR #69032:
URL: https://github.com/apache/doris/pull/69032#discussion_r4381111390
##########
be/src/exprs/function/geo/functions_geo.cpp:
##########
@@ -18,32 +18,1301 @@
#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 std::unique_ptr<GeoShape> decode_geography_calculation_shape(StringRef
value);
+
+static Status decode_calculation_shape(StringRef value, const DataTypePtr&
type,
+ std::unique_ptr<GeoShape>* shape,
+ const char* function_name) {
+ *shape = is_spatial_type(type) ? decode_geography_calculation_shape(value)
+ : decode_geo_shape(value, type);
+ // Typed admission rejects unsupported and EMPTY calculation shapes before
distance
+ // helpers can expose infinity/max sentinels. Raw WKB storage/accessors
remain unchanged.
+ if (is_spatial_type(type) && *shape == nullptr) {
+ 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);
+ }
+};
+
+// Raw identity must survive ToPoint rounding: 58 and 58.00000000000001 degrees
+// produce the same S2Point, but only one is an endpoint of the original edge.
+struct GeographyVertex {
+ GeoCoordinate coordinate;
+ S2Point point;
+
+ GeographyVertex() : GeographyVertex({.x = 0, .y = 0}) {}
+ explicit GeographyVertex(GeoCoordinate value) : coordinate(value) {
+ if (std::abs(coordinate.y) == 90) {
+ coordinate.x = 0;
+ point = {0, 0, coordinate.y > 0 ? 1.0 : -1.0};
+ } else {
+ if (coordinate.x == 180) {
+ coordinate.x = -180;
+ }
+ point = S2LatLng::FromDegrees(coordinate.y,
coordinate.x).ToPoint();
+ }
+ }
+ bool operator==(const GeographyVertex& other) const {
+ return coordinate.x == other.coordinate.x && coordinate.y ==
other.coordinate.y;
+ }
+ bool operator<(const GeographyVertex& other) const {
+ return std::pair {coordinate.x, coordinate.y} <
+ std::pair {other.coordinate.x, other.coordinate.y};
+ }
+};
+
+static bool geography_antipodal(const GeographyVertex& a, const
GeographyVertex& b) {
+ return a.coordinate.y == -b.coordinate.y &&
+ (std::abs(a.coordinate.y) == 90 ||
+ geography_degree_sum(std::max(a.coordinate.x, b.coordinate.x),
+ -std::min(a.coordinate.x, b.coordinate.x)) ==
+ std::pair {180.0, 0.0});
+}
+
+struct GeographyCoordinates {
+ GeoShapeType type = GeoShapeType::GEO_SHAPE_ANY;
+ DorisVector<DorisVector<GeographyVertex>> sequences;
+ GeoParseStatus status;
+
+ explicit GeographyCoordinates(StringRef wkb) {
+ GeoCoordinateListList raw;
+ status = WkbParse::calculation_coordinates(wkb.data, wkb.size, &raw,
&type);
+ if (status != GeoParseStatus::GEO_PARSE_OK) {
+ return;
+ }
+ for (const auto& sequence : raw.list) {
+ auto& vertices = sequences.emplace_back();
+ for (const auto& coordinate : sequence->list) {
+ vertices.emplace_back(coordinate);
+ }
+ }
+ }
+};
+
+static bool geography_stable_points(const DorisVector<GeographyVertex>&
sequence,
+ std::vector<S2Point>* points) {
+ for (size_t i = 0; i + 1 < sequence.size(); ++i) {
+ const auto& a = sequence[i];
+ const auto& b = sequence[i + 1];
+ // (0,0)->(180,0) has no unique arc and would normalize a zero raw
normal.
+ if (geography_antipodal(a, b)) {
+ return false;
+ }
+ points->push_back(a.point);
+ // Stabilize before S2 admission: rounded (30,58)->(-150,-58-e) can be
+ // exactly antipodal despite its unique raw arc passing through the
south pole.
+ if (a.point.DotProd(b.point) < 0) {
+ points->push_back(GeographyRawArc(a.coordinate,
b.coordinate).midpoint);
+ }
+ }
+ // A line's final vertex is real, including its closing edge; only rings
remove it.
+ points->push_back(sequence.back().point);
+ points->erase(std::ranges::unique(*points).begin(), points->end());
+ return true;
+}
+
+static std::unique_ptr<GeoShape> decode_geography_line(const
GeographyCoordinates& raw,
+ StringRef value) {
+ std::vector<S2Point> points;
+ if (!geography_stable_points(raw.sequences.front(), &points)) {
+ return nullptr;
+ }
+ if (points.size() < 2) {
+ // Canonical pole aliases such as (0,90),(90,90) collapse to one
calculation
+ // point. Preserve their prior admission; the old parser still rejects
a
+ // truly duplicate-only line such as (1,1),(1,1).
+ GeoParseStatus status;
+ return GeoShape::from_wkb_bytes(value.data, value.size, status);
+ }
+ auto polyline = std::make_unique<S2Polyline>(points, S2Debug::DISABLE);
+ if (!polyline->IsValid()) {
+ return nullptr;
+ }
+ return GeoLine::create_unique(std::move(polyline));
+}
+
+static std::unique_ptr<GeoShape> decode_geography_calculation_shape(StringRef
value) {
+ GeographyCoordinates raw(value);
+ if (raw.status != GeoParseStatus::GEO_PARSE_OK) {
+ return nullptr;
+ }
+ if (raw.type == GeoShapeType::GEO_SHAPE_POINT) {
+ GeoParseStatus status;
+ return GeoShape::from_wkb_bytes(value.data, value.size, status);
+ }
+ if (raw.type == GeoShapeType::GEO_SHAPE_LINE_STRING) {
+ return decode_geography_line(raw, value);
+ }
+ std::vector<std::unique_ptr<S2Loop>> loops;
+ std::vector<const S2Loop*> input_loops;
+ for (const auto& ring : raw.sequences) {
+ std::vector<S2Point> points;
+ // Validate holes against stable raw arcs, as with typed line
admission.
+ if (!geography_stable_points(ring, &points)) {
+ return nullptr;
+ }
+ // Validated WKB rings are closed; only polygon admission removes that
endpoint.
+ DCHECK(points.front() == points.back());
+ points.pop_back();
+ if (points.size() < 3) {
+ return nullptr;
+ }
+ auto loop = std::make_unique<S2Loop>(points, S2Debug::DISABLE);
+ if (!loop->IsValid()) {
+ return nullptr;
+ }
+ loop->Normalize();
+ if (!loops.empty() && !loops.front()->Contains(*loop)) {
+ return nullptr;
+ }
+ input_loops.push_back(loop.get());
+ loops.push_back(std::move(loop));
+ }
+ auto polygon = std::make_unique<S2Polygon>(std::move(loops),
S2Debug::DISABLE);
+ if (!polygon->IsValid()) {
+ return nullptr;
+ }
+ std::vector<int> depths;
+ for (const auto* loop : input_loops) {
+ depths.push_back(loop->depth());
+ }
+ return GeoPolygon::create_unique(std::move(polygon), std::move(depths));
+}
+
+#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:
+ GeographyCoordinates left;
+ GeographyCoordinates right;
+
+ GeographyEdgeNodes(StringRef left_wkb, StringRef right_wkb) :
left(left_wkb), right(right_wkb) {
+ DCHECK_EQ(left.status, GeoParseStatus::GEO_PARSE_OK);
+ DCHECK_EQ(right.status, GeoParseStatus::GEO_PARSE_OK);
+ DorisMap<double, DorisVector<GeographyVertex>> longitude_vertices;
+ std::set<GeographyVertex, std::less<>,
CustomStdAllocator<GeographyVertex>> coordinates;
+ for (const auto* shape : {&left, &right}) {
+ for (const auto& sequence : shape->sequences) {
+ for (const auto& point : sequence) {
+ if (!coordinates.insert(point).second) {
+ continue;
+ }
+ const auto& coordinate = point.coordinate;
+ longitude_vertices[coordinate.x].push_back(point);
+ if (std::abs(coordinate.y) != 90) {
+ _meridians[meridian(coordinate.x)][coordinate.x <
0].push_back(
+ {coordinate.y, point});
+ }
+ }
+ }
+ }
+ // Keep dense longitude groups separately so a circle can skip its
known
+ // endpoint longitudes before enumerating nearby points. At sqrt(N)
points
+ // per group there are at most sqrt(N) extra indexes, even with
unrelated
+ // vertices padding a near-antipodal spoke input.
+ const auto dense_size =
+ std::max(size_t {2},
static_cast<size_t>(std::ceil(std::sqrt(coordinates.size()))));
+ for (const auto& [longitude, points] : longitude_vertices) {
+ auto& index =
+ points.size() >= dense_size ? _dense_longitudes[longitude]
: _vertex_index;
+ for (const auto& point : points) {
+ index.Add(point.point, point);
+ }
+ }
+ for (auto& entry : _meridians) {
+ for (auto& side : entry.second) {
+ side.push_back({-90, GeographyVertex({.x = 0, .y = -90})});
+ side.push_back({90, GeographyVertex({.x = 0, .y = 90})});
+ std::ranges::sort(side, {}, &Vertex::latitude);
+ }
+ }
+ }
+
+ bool contains_point(const GeographyVertex& start, const GeographyVertex&
end,
+ const GeographyVertex& point) const {
+ // Preserve exact endpoints, including date-line and pole aliases.
+ if (point == start || point == end) {
+ return true;
+ }
+ if (start == end) {
+ return false;
+ }
+ const auto& a = start.coordinate;
+ const auto& b = end.coordinate;
+ const auto& c = point.coordinate;
+ const auto key = meridian(std::abs(a.y) == 90 ? b.x : a.x);
+ bool member = false;
+ if (meridian_edge(start, end, [&](const auto& vertices, double from,
double to) {
+ // Opposite-longitude edges can pass through a pole; the pole
belongs
+ // to both legs, whereas other positions must be on the
selected side.
+ const bool same_side =
+ std::abs(c.y) == 90 ||
+ (meridian(c.x) == key && &_meridians.at(key)[c.x < 0]
== &vertices);
+ member |= same_side && c.y >= std::min(from, to) && c.y <=
std::max(from, to);
+ })) {
+ return member;
+ }
+ const bool wraps = geography_degree_sum(std::max(a.x, b.x),
-std::min(a.x, b.x)) >
+ std::pair {180.0, 0.0};
+ const bool within = wraps ? (c.x < std::min(a.x, b.x) || c.x >
std::max(a.x, b.x))
+ : (c.x > std::min(a.x, b.x) && c.x <
std::max(a.x, b.x));
+ // Near-antipodal endpoint conversion can rotate the S2 plane
substantially.
+ // POINT(45 0) is on (-44.99999999999997 -45,134.99999999999997 45);
+ // decide from raw degrees rather than accepting a larger distance
tolerance.
+ return within && geography_common_great_circle(a, b, c);
+ }
+
+ void append_edge(const GeographyVertex& start, const GeographyVertex& end,
+ DorisVector<GeographyVertex>* output) const {
+ output->push_back(start);
+ if (start == end) {
+ return;
+ }
+ if (meridian_edge(start, end, [&](const auto& vertices, double from,
double to) {
+ append_range(vertices, from, to, output);
+ if (std::abs(to) == 90) {
+ const GeographyVertex pole({.x = 0, .y = to});
+ if (end != pole) {
+ output->push_back(pole);
+ }
+ }
+ })) {
+ return;
+ }
+ const auto& a = start.coordinate;
+ const auto& b = end.coordinate;
+ const bool wraps = geography_degree_sum(std::max(a.x, b.x),
-std::min(a.x, b.x)) >
+ std::pair {180.0, 0.0};
+ const bool east = (a.x < b.x) != wraps;
+ // A near-antipodal rounded edge can miss true raw-coordinate nodes.
Two
+ // stable half-arcs recall those nodes without a hemispherical search.
+ const auto midpoint = GeographyRawArc(a, b).midpoint;
+ DorisVector<GeographyVertex> on_edge;
+ for (const auto& endpoint : {start, end}) {
+ visit_edge_candidates(
+ endpoint.point, midpoint, [](double) { return false; },
+ [&](const GeographyVertex& point) {
+ const auto& c = point.coordinate;
+ // Longitude is strictly monotone along every
non-meridian minor arc.
+ // Check the original interval (including date-line
wrap), not a distance
+ // tolerance that could append a point just beyond an
endpoint.
+ const bool within =
+ wraps ? (c.x < std::min(a.x, b.x) || c.x >
std::max(a.x, b.x))
+ : (c.x > std::min(a.x, b.x) && c.x <
std::max(a.x, b.x));
+ if (within && geography_common_great_circle(a, b, c)) {
+ on_edge.push_back(point);
+ }
+ });
+ }
+ std::ranges::sort(on_edge, [&](const GeographyVertex& x, const
GeographyVertex& y) {
+ const double x_lon = x.coordinate.x;
+ const double y_lon = y.coordinate.x;
+ if ((x_lon < a.x) != (y_lon < a.x)) {
+ return east ? x_lon >= a.x : x_lon < a.x;
+ }
+ return east ? x_lon < y_lon : x_lon > y_lon;
+ });
+ on_edge.erase(std::ranges::unique(on_edge).begin(), on_edge.end());
+ output->insert(output->end(), on_edge.begin(), on_edge.end());
+ }
+
+ // Calculation-only points never enter the raw-degree identity/noding map.
+ std::optional<S2Point> arc_midpoint(const GeographyVertex& a, const
GeographyVertex& b) const {
+ if (a.point.DotProd(b.point) < 0) {
+ return GeographyRawArc(a.coordinate, b.coordinate).midpoint;
+ }
+ return std::nullopt;
+ }
+
+ void add_line(const GeographyCoordinates& shape, MutableS2ShapeIndex*
index) const {
+ using Edge = std::pair<GeographyVertex, GeographyVertex>;
+ using EdgeSet = std::set<Edge, std::less<>, CustomStdAllocator<Edge>>;
+ using Interval = std::pair<double, double>;
+ DorisMap<const DorisVector<Vertex>*, DorisVector<Interval>> intervals;
+ EdgeSet raw_edges;
+ const auto& line = shape.sequences.front();
+ const auto& first = line.front();
+ const auto& last = line.back();
+ for (size_t i = 0; i + 1 < line.size(); ++i) {
+ const auto& a = line[i];
+ const auto& b = line[i + 1];
+ if (a == b) {
+ continue;
+ }
+ if (!meridian_edge(a, b, [&](const auto& vertices, double from,
double to) {
+ intervals[&vertices].emplace_back(std::min(from, to),
std::max(from, to));
+ })) {
+ raw_edges.emplace(std::min(a, b), std::max(a, b));
+ }
+ }
+ EdgeSet edges;
+ auto add_edge = [&](const GeographyVertex& a, const GeographyVertex&
b) {
+ edges.emplace(std::min(a, b), std::max(a, b));
+ };
+ // Merge latitude coverage before splitting: N repeated or overlapping
edges
+ // and M collinear vertices produce O(N+M) scratch, never an N*M
vertex walk.
+ for (auto& [vertices, ranges] : intervals) {
+ std::ranges::sort(ranges);
+ auto emit = [&](double from, double to) {
+ auto begin = std::ranges::lower_bound(*vertices, from, {},
&Vertex::latitude);
+ const auto end = std::ranges::upper_bound(*vertices, to, {},
&Vertex::latitude);
+ for (auto next = begin + 1; next != end; ++begin, ++next) {
+ add_edge(begin->point, next->point);
+ }
+ };
+ double from = ranges.front().first;
+ double to = ranges.front().second;
+ for (const auto& range : ranges) {
+ if (range.first > to) {
+ emit(from, to);
+ from = range.first;
+ }
+ to = std::max(to, range.second);
+ }
+ emit(from, to);
+ }
+ add_circle_edges(raw_edges, add_edge);
+ auto add = [&](const S2Point& a, const S2Point& b) {
+ index->Add(std::make_unique<GeographyLineShape>(a, b));
+ };
+ if (edges.empty()) {
+ // A raw pole-alias line may collapse to one canonical calculation
position.
+ add(first.point, first.point);
+ }
+ std::set<S2Point, std::less<>, CustomStdAllocator<S2Point>>
interior_vertices;
+ for (const auto& [a, b] : edges) {
+ const auto midpoint = arc_midpoint(a, b);
+ if (midpoint) {
+ add(a.point, *midpoint);
+ add(*midpoint, b.point);
+ // The synthetic split is interior, not a new OPEN line
boundary:
+ // a crossing at this point must not become ST_Touches=true.
+ interior_vertices.insert(*midpoint);
+ } else {
+ add(a.point, b.point);
Review Comment:
[P2] Preserve open endpoints when atomic edges round to zero length
For default GEOGRAPHY values imported as WKB, let A=`LINESTRING(0 58,0 59)`,
A2=`LINESTRING(0 58,0 58.00000000000001,0 59)` and B=`LINESTRING(-1 57,0 58,1
57)`. These are valid nondegenerate XY lines; A and A2 describe the same
meridian arc and meet B only at their original boundary endpoint, so ST_Touches
must be true for both, in either operand order.
Pinned S2's [degree
conversion](https://github.com/google/s2geometry/blob/v0.10.0/src/s2/s1angle.h#L290-L292)
maps both latitudes 58 and 58.00000000000001 to the same radians
`0x1.03258031ade7cp+0`, hence identical S2Points. Admission removes that
adjacent rounded duplicate, but the raw meridian sweep retains it and this
unconditional emission creates an AA edge at the original endpoint. B also
receives an AA marker for its genuine interior vertex there. S2 [retains
matching degenerate polylines in OPEN
intersection](https://github.com/google/s2geometry/blob/v0.10.0/src/s2/s2boolean_operation.cc#L1397-L1406),
while an ordinary OPEN arc excludes its endpoint. Consequently the
OPEN-intersection test changes from false to true and ST_Touches changes from
true to false solely under subdivision.
This is a downstream nonpoint boundary failure, distinct from repaired
raw-map identity/antipodal issue 4379422808 and non-meridian equality issue
4374330645. Preserve the original boundary role when handling rounded-zero
atomic edges, while retaining genuine interior markers; cover
subdivision/reversal and both operand orders. Evidence is static
source/arithmetic analysis, not a Doris runtime execution.
--
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]
