OIiveirra commented on code in PR #69032: URL: https://github.com/apache/doris/pull/69032#discussion_r4404258153
########## 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: Fixed in 130b134012e. Keep native wide arithmetic where sufficient and use Boost software quad on binary64-long-double targets, without weakening the precision/exponent contract. Native/53-bit-ABI/GCC15 compiler probes and forced software arithmetic tests passed; a full Apple SDK build remains for platform CI. ########## 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: Fixed in 130b134012e. Tiny Polygon edges can use complete indexed raw-coordinate bands; larger or unproven bands retain the full S2 fallback. Exact membership and directed rings/holes are unchanged, with no truncation or new input limits. Dense/axis-swapped work-count tests, complete-fallback controls, 152 BE tests and three targeted native regression suites passed. -- 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]
