github-actions[bot] commented on code in PR #69032: URL: https://github.com/apache/doris/pull/69032#discussion_r4402847165
########## be/src/exprs/function/geo/geography_calculation.cpp: ########## @@ -0,0 +1,1402 @@ +// 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 "exprs/function/geo/geography_calculation.h" + +#include <glog/logging.h> +#include <s2/mutable_s2shape_index.h> +#include <s2/s2boolean_operation.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> +// cpp_int needs no quadmath wrappers. Boost's float128 numeric_limits specialization +// conflicts with libstdc++ 15's native specialization in the GCC performance build. +#if defined(__GNUC__) && !defined(__clang__) && __GNUC__ >= 15 +#define BOOST_CSTDFLOAT_NO_LIBQUADMATH_SUPPORT +#endif +#include <boost/multiprecision/cpp_int.hpp> +#if defined(__GNUC__) && !defined(__clang__) && __GNUC__ >= 15 +#undef BOOST_CSTDFLOAT_NO_LIBQUADMATH_SUPPORT +#endif +#include <cmath> +#include <cstdint> +#include <limits> +#include <map> +#include <numbers> +#include <optional> +#include <set> +#include <utility> + +#include "core/assert_cast.h" +#include "core/custom_allocator.h" +#include "exprs/function/geo/geo_types.h" +#include "exprs/function/geo/wkb_parse.h" + +namespace doris { + +// 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 && Review Comment: [P1] Preserve Apple ARM64 builds when selecting wide arithmetic Apple's [ARM64 ABI](https://developer.apple.com/documentation/xcode/writing-arm64-code-for-apple-platforms?changes=_1) defines `long double` identically to `double`, so its `digits` is 53 and `min_exponent` is -1021. This unconditional assertion therefore fails on Apple Silicon. The repository explicitly supports Darwin/ARM in `be/CMakeLists.txt`, and `be/src/exprs/CMakeLists.txt` includes this new translation unit unconditionally, so BE and BE-unit-test builds fail even when no Geography function is used. Use a portable wider representation or equivalent arithmetic that retains the required residual precision/exponent range on that target; removing the assertion alone would lose the guarantees described here. This is a source/ABI proof; no build was run. ########## be/src/exprs/function/geo/geography_calculation.cpp: ########## @@ -0,0 +1,1402 @@ +// 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 "exprs/function/geo/geography_calculation.h" + +#include <glog/logging.h> +#include <s2/mutable_s2shape_index.h> +#include <s2/s2boolean_operation.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> +// cpp_int needs no quadmath wrappers. Boost's float128 numeric_limits specialization +// conflicts with libstdc++ 15's native specialization in the GCC performance build. +#if defined(__GNUC__) && !defined(__clang__) && __GNUC__ >= 15 +#define BOOST_CSTDFLOAT_NO_LIBQUADMATH_SUPPORT +#endif +#include <boost/multiprecision/cpp_int.hpp> +#if defined(__GNUC__) && !defined(__clang__) && __GNUC__ >= 15 +#undef BOOST_CSTDFLOAT_NO_LIBQUADMATH_SUPPORT +#endif +#include <cmath> +#include <cstdint> +#include <limits> +#include <map> +#include <numbers> +#include <optional> +#include <set> +#include <utility> + +#include "core/assert_cast.h" +#include "core/custom_allocator.h" +#include "exprs/function/geo/geo_types.h" +#include "exprs/function/geo/wkb_parse.h" + +namespace doris { + +// 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)); +} + +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. +// Keep this exact cancellation proof together during the behavior-preserving module extraction. +// NOLINTNEXTLINE(readability-function-size) +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: + using LineEdge = std::pair<GeographyVertex, GeographyVertex>; + using LineEdges = std::set<LineEdge, std::less<>, CustomStdAllocator<LineEdge>>; + using PolygonEdges = std::array<LineEdges, 2>; + + struct RawEdges { + LineEdges line; + PolygonEdges polygon; + + auto vertices(const GeographyCoordinates& raw) const { + std::set<GeographyVertex, std::less<>, CustomStdAllocator<GeographyVertex>> result; + // A pole-alias line can have no edges but still denotes one canonical position. + result.insert(raw.sequences.front().front()); + const auto append = [&](const LineEdges& edges) { + for (const auto& [a, b] : edges) { + result.insert(a); + result.insert(b); + } + }; + append(line); + for (const auto& ring_edges : polygon) { + append(ring_edges); + } + return result; + } + }; + + 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); + } + + bool line_contains_point(const GeographyCoordinates& shape, + const GeographyVertex& point) const { + const auto& line = shape.sequences.front(); + for (size_t i = 0; i + 1 < line.size(); ++i) { + if (contains_point(line[i], line[i + 1], point)) { + return true; + } + } + return false; + } + + bool line_boundary_touches(const GeographyCoordinates& line, + const GeographyCoordinates& other) const { + const auto& vertices = line.sequences.front(); + // Closed lines have no boundary. An open endpoint must belong to the other raw line: + // a rounded coincidence across the gap 58 -> 58+ulp is not a contact. + return vertices.front() != vertices.back() && + (line_contains_point(other, vertices.front()) || + line_contains_point(other, vertices.back())); + } + + 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)) { Review Comment: [P2] Bound exact noding work for dense Polygon vertices Polygon rings call `append_edge` for every original edge, bypassing the LineString batching. A concrete default-Geography WKB ring uses even N=10,000, r=1e-14 and vertices `((-r if i%2==0 else r)+0.1*r*i/(N-1), -r+2*r*i/(N-1))` for i=0..N-1, then `(2*r,2*r),(2*r,-2*r),first`. It is a simple, nondegenerate ring: the zigzag rises monotonically and the closing edges stay outside it. All converted vertices lie within 1e-15 radians of every edge endpoint, below the fixed 64*epsilon query radius; the [pinned query returns all candidates within that radius](https://github.com/google/s2geometry/blob/v0.10.0/src/s2/s2closest_point_query.h#L27-L68). Both queries per edge therefore enumerate the whole cluster. This strict-longitude gate causes N*N-2*N+2 = **99,980,002** calls to the 48-term bigint/map proof for only **160,061 bytes** before `ST_Area_Square_Meters` can return. Final noding output remains linear. Bound or batch Polygon candidate/proof work and add increasing-size work-count coverage. This is the Polygon per-edge consumer, separate from existing 4375245664/4376498681's repaired LineString circle paths. Counts follow from source/arithmetic modeling; no 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]
