|
25 | 25 | #include <s2/s2cap.h> |
26 | 26 | #include <s2/s2earth.h> |
27 | 27 | #include <s2/s2edge_crosser.h> |
| 28 | +#include <s2/s2edge_distances.h> |
28 | 29 | #include <s2/s2latlng.h> |
29 | 30 | #include <s2/s2loop.h> |
30 | 31 | #include <s2/s2point.h> |
|
38 | 39 | // IWYU pragma: no_include <bits/std_abs.h> |
39 | 40 | #include <cmath> |
40 | 41 | #include <iomanip> |
| 42 | +#include <limits> |
41 | 43 | #include <sstream> |
42 | 44 | #include <utility> |
43 | 45 | #include <vector> |
@@ -1630,5 +1632,374 @@ std::string GeoShape::as_binary(GeoShape* rhs) { |
1630 | 1632 | return res; |
1631 | 1633 | } |
1632 | 1634 |
|
| 1635 | +// Helper function to compute distance from a point to a line segment |
| 1636 | +double distance_point_to_segment(const S2Point& point, const S2Point& line_start, |
| 1637 | + const S2Point& line_end) { |
| 1638 | + S1Angle angle = S2::GetDistance(point, line_start, line_end); |
| 1639 | + return S2Earth::ToMeters(angle); |
| 1640 | +} |
| 1641 | + |
| 1642 | +// Helper function to compute distance from a point to a polyline |
| 1643 | +double distance_point_to_polyline(const S2Point& point, const S2Polyline* polyline) { |
| 1644 | + if (polyline->num_vertices() == 0) { |
| 1645 | + return std::numeric_limits<double>::max(); |
| 1646 | + } |
| 1647 | + if (polyline->num_vertices() == 1) { |
| 1648 | + return S2Earth::GetDistanceMeters(point, polyline->vertex(0)); |
| 1649 | + } |
| 1650 | + |
| 1651 | + S1Angle min_angle = S1Angle::Infinity(); |
| 1652 | + for (int i = 0; i < polyline->num_vertices() - 1; ++i) { |
| 1653 | + const S2Point& p1 = polyline->vertex(i); |
| 1654 | + const S2Point& p2 = polyline->vertex(i + 1); |
| 1655 | + |
| 1656 | + S1Angle dist = S2::GetDistance(point, p1, p2); |
| 1657 | + if (dist < min_angle) { |
| 1658 | + min_angle = dist; |
| 1659 | + } |
| 1660 | + } |
| 1661 | + |
| 1662 | + return S2Earth::ToMeters(min_angle); |
| 1663 | +} |
| 1664 | + |
| 1665 | +// Helper function to compute distance from a point to a polygon |
| 1666 | +double distance_point_to_polygon(const S2Point& point, const S2Polygon* polygon) { |
| 1667 | + // Check if point is inside polygon |
| 1668 | + if (polygon->Contains(point)) { |
| 1669 | + return 0.0; |
| 1670 | + } |
| 1671 | + |
| 1672 | + // Find minimum distance to polygon boundary |
| 1673 | + S1Angle min_angle = S1Angle::Infinity(); |
| 1674 | + |
| 1675 | + for (int i = 0; i < polygon->num_loops(); ++i) { |
| 1676 | + const S2Loop* loop = polygon->loop(i); |
| 1677 | + |
| 1678 | + for (int j = 0; j < loop->num_vertices(); ++j) { |
| 1679 | + const S2Point& p1 = loop->vertex(j); |
| 1680 | + const S2Point& p2 = loop->vertex((j + 1) % loop->num_vertices()); |
| 1681 | + |
| 1682 | + S1Angle dist = S2::GetDistance(point, p1, p2); |
| 1683 | + if (dist < min_angle) { |
| 1684 | + min_angle = dist; |
| 1685 | + } |
| 1686 | + } |
| 1687 | + } |
| 1688 | + |
| 1689 | + return S2Earth::ToMeters(min_angle); |
| 1690 | +} |
| 1691 | + |
| 1692 | +double GeoPoint::Length() const { |
| 1693 | + // Point has no length |
| 1694 | + return 0.0; |
| 1695 | +} |
| 1696 | + |
| 1697 | +double GeoLine::Length() const { |
| 1698 | + // GeoLine is always valid with at least 2 vertices (guaranteed by constructor) |
| 1699 | + double total_length = 0.0; |
| 1700 | + for (int i = 0; i < _polyline->num_vertices() - 1; ++i) { |
| 1701 | + const S2Point& p1 = _polyline->vertex(i); |
| 1702 | + const S2Point& p2 = _polyline->vertex(i + 1); |
| 1703 | + |
| 1704 | + S2LatLng lat_lng1(p1); |
| 1705 | + S2LatLng lat_lng2(p2); |
| 1706 | + |
| 1707 | + // Calculate distance in meters using S2Earth |
| 1708 | + double distance = S2Earth::GetDistanceMeters(lat_lng1, lat_lng2); |
| 1709 | + total_length += distance; |
| 1710 | + } |
| 1711 | + |
| 1712 | + return total_length; |
| 1713 | +} |
| 1714 | + |
| 1715 | +double GeoPolygon::Length() const { |
| 1716 | + // GeoPolygon is always valid with at least one loop (guaranteed by constructor) |
| 1717 | + double perimeter = 0.0; |
| 1718 | + |
| 1719 | + for (int loop_idx = 0; loop_idx < _polygon->num_loops(); ++loop_idx) { |
| 1720 | + const S2Loop* loop = _polygon->loop(loop_idx); |
| 1721 | + for (int i = 0; i < loop->num_vertices(); ++i) { |
| 1722 | + const S2Point& p1 = loop->vertex(i); |
| 1723 | + const S2Point& p2 = loop->vertex((i + 1) % loop->num_vertices()); |
| 1724 | + |
| 1725 | + S2LatLng lat_lng1(p1); |
| 1726 | + S2LatLng lat_lng2(p2); |
| 1727 | + |
| 1728 | + // Calculate distance in meters using S2Earth |
| 1729 | + double distance = S2Earth::GetDistanceMeters(lat_lng1, lat_lng2); |
| 1730 | + perimeter += distance; |
| 1731 | + } |
| 1732 | + } |
| 1733 | + |
| 1734 | + return perimeter; |
| 1735 | +} |
| 1736 | + |
| 1737 | +double GeoMultiPolygon::Length() const { |
| 1738 | + double total_length = 0.0; |
| 1739 | + |
| 1740 | + // Calculate the perimeter of all polygons |
| 1741 | + for (const auto& polygon : _polygons) { |
| 1742 | + total_length += polygon->Length(); |
| 1743 | + } |
| 1744 | + |
| 1745 | + return total_length; |
| 1746 | +} |
| 1747 | + |
| 1748 | +double GeoCircle::Length() const { |
| 1749 | + // GeoCircle is always valid (guaranteed by constructor) |
| 1750 | + // Get the radius in meters |
| 1751 | + double radius_meters = S2Earth::ToMeters(_cap->radius()); |
| 1752 | + |
| 1753 | + // Calculate circumference: 2 * pi * r |
| 1754 | + return 2.0 * M_PI * radius_meters; |
| 1755 | +} |
| 1756 | + |
| 1757 | +double GeoPoint::Distance(const GeoShape* rhs) const { |
| 1758 | + // rhs is guaranteed to be valid by StDistance (functions_geo.cpp) |
| 1759 | + switch (rhs->type()) { |
| 1760 | + case GEO_SHAPE_POINT: { |
| 1761 | + const GeoPoint* point = static_cast<const GeoPoint*>(rhs); |
| 1762 | + S2LatLng this_ll = S2LatLng(*_point); |
| 1763 | + S2LatLng other_ll = S2LatLng(*point->point()); |
| 1764 | + return S2Earth::GetDistanceMeters(this_ll, other_ll); |
| 1765 | + } |
| 1766 | + case GEO_SHAPE_LINE_STRING: { |
| 1767 | + const GeoLine* line = static_cast<const GeoLine*>(rhs); |
| 1768 | + return distance_point_to_polyline(*_point, line->polyline()); |
| 1769 | + } |
| 1770 | + case GEO_SHAPE_POLYGON: { |
| 1771 | + const GeoPolygon* polygon = static_cast<const GeoPolygon*>(rhs); |
| 1772 | + return distance_point_to_polygon(*_point, polygon->polygon()); |
| 1773 | + } |
| 1774 | + case GEO_SHAPE_CIRCLE: { |
| 1775 | + const GeoCircle* circle = static_cast<const GeoCircle*>(rhs); |
| 1776 | + S2LatLng this_ll = S2LatLng(*_point); |
| 1777 | + S2LatLng center_ll = S2LatLng(circle->circle()->center()); |
| 1778 | + double dist_to_center = S2Earth::GetDistanceMeters(this_ll, center_ll); |
| 1779 | + double circle_radius = S2Earth::ToMeters(circle->circle()->radius()); |
| 1780 | + |
| 1781 | + // Distance from point to circle is distance to center minus radius |
| 1782 | + if (dist_to_center <= circle_radius + TOLERANCE) { |
| 1783 | + return 0.0; // Point is inside or on the circle boundary |
| 1784 | + } |
| 1785 | + return dist_to_center - circle_radius; |
| 1786 | + } |
| 1787 | + case GEO_SHAPE_MULTI_POLYGON: { |
| 1788 | + return rhs->Distance(this); // Delegate to MultiPolygon's implementation |
| 1789 | + } |
| 1790 | + default: |
| 1791 | + return -1.0; |
| 1792 | + } |
| 1793 | +} |
| 1794 | + |
| 1795 | +double GeoLine::Distance(const GeoShape* rhs) const { |
| 1796 | + // rhs is guaranteed to be valid by StDistance (functions_geo.cpp) |
| 1797 | + switch (rhs->type()) { |
| 1798 | + case GEO_SHAPE_POINT: { |
| 1799 | + return rhs->Distance(this); // Reuse Point's Distance implementation |
| 1800 | + } |
| 1801 | + case GEO_SHAPE_LINE_STRING: { |
| 1802 | + const GeoLine* other_line = static_cast<const GeoLine*>(rhs); |
| 1803 | + if (_polyline->Intersects(*other_line->polyline())) { |
| 1804 | + return 0.0; |
| 1805 | + } |
| 1806 | + double min_distance = std::numeric_limits<double>::max(); |
| 1807 | + |
| 1808 | + // Check distance from each vertex of this line to other line |
| 1809 | + for (int i = 0; i < _polyline->num_vertices(); ++i) { |
| 1810 | + double dist = distance_point_to_polyline(_polyline->vertex(i), other_line->polyline()); |
| 1811 | + min_distance = std::min(min_distance, dist); |
| 1812 | + } |
| 1813 | + |
| 1814 | + // Check distance from each vertex of other line to this line |
| 1815 | + for (int i = 0; i < other_line->polyline()->num_vertices(); ++i) { |
| 1816 | + double dist = |
| 1817 | + distance_point_to_polyline(other_line->polyline()->vertex(i), _polyline.get()); |
| 1818 | + min_distance = std::min(min_distance, dist); |
| 1819 | + } |
| 1820 | + |
| 1821 | + // Handle touching case: if min_distance is within tolerance, lines are touching |
| 1822 | + if (min_distance <= TOLERANCE) { |
| 1823 | + return 0.0; |
| 1824 | + } |
| 1825 | + return min_distance; |
| 1826 | + } |
| 1827 | + case GEO_SHAPE_POLYGON: { |
| 1828 | + return rhs->Distance(this); // Delegate to Polygon's implementation |
| 1829 | + } |
| 1830 | + case GEO_SHAPE_CIRCLE: { |
| 1831 | + return rhs->Distance(this); // Delegate to Circle's implementation |
| 1832 | + } |
| 1833 | + case GEO_SHAPE_MULTI_POLYGON: { |
| 1834 | + return rhs->Distance(this); // Delegate to MultiPolygon's implementation |
| 1835 | + } |
| 1836 | + default: |
| 1837 | + return -1.0; |
| 1838 | + } |
| 1839 | +} |
| 1840 | + |
| 1841 | +double GeoPolygon::Distance(const GeoShape* rhs) const { |
| 1842 | + // rhs is guaranteed to be valid by StDistance (functions_geo.cpp) |
| 1843 | + switch (rhs->type()) { |
| 1844 | + case GEO_SHAPE_POINT: { |
| 1845 | + return rhs->Distance(this); // Reuse Point's Distance implementation |
| 1846 | + } |
| 1847 | + case GEO_SHAPE_LINE_STRING: { |
| 1848 | + const GeoLine* line = static_cast<const GeoLine*>(rhs); |
| 1849 | + if (_polygon->Intersects(*line->polyline())) { |
| 1850 | + return 0.0; |
| 1851 | + } |
| 1852 | + double min_distance = std::numeric_limits<double>::max(); |
| 1853 | + |
| 1854 | + // Check distance from each vertex of line to polygon |
| 1855 | + for (int i = 0; i < line->polyline()->num_vertices(); ++i) { |
| 1856 | + double dist = distance_point_to_polygon(line->polyline()->vertex(i), _polygon.get()); |
| 1857 | + min_distance = std::min(min_distance, dist); |
| 1858 | + } |
| 1859 | + |
| 1860 | + // Check distance from each polygon vertex to line |
| 1861 | + for (int i = 0; i < _polygon->num_loops(); ++i) { |
| 1862 | + const S2Loop* loop = _polygon->loop(i); |
| 1863 | + for (int j = 0; j < loop->num_vertices(); ++j) { |
| 1864 | + double dist = distance_point_to_polyline(loop->vertex(j), line->polyline()); |
| 1865 | + min_distance = std::min(min_distance, dist); |
| 1866 | + } |
| 1867 | + } |
| 1868 | + |
| 1869 | + // Handle touching case: if min_distance is within tolerance, they are touching |
| 1870 | + if (min_distance <= TOLERANCE) { |
| 1871 | + return 0.0; |
| 1872 | + } |
| 1873 | + return min_distance; |
| 1874 | + } |
| 1875 | + case GEO_SHAPE_POLYGON: { |
| 1876 | + const GeoPolygon* other = static_cast<const GeoPolygon*>(rhs); |
| 1877 | + if (_polygon->Intersects(*other->polygon())) { |
| 1878 | + return 0.0; |
| 1879 | + } |
| 1880 | + double min_distance = std::numeric_limits<double>::max(); |
| 1881 | + |
| 1882 | + // Check distance from each vertex of this polygon to other polygon |
| 1883 | + for (int i = 0; i < _polygon->num_loops(); ++i) { |
| 1884 | + const S2Loop* loop = _polygon->loop(i); |
| 1885 | + for (int j = 0; j < loop->num_vertices(); ++j) { |
| 1886 | + double dist = distance_point_to_polygon(loop->vertex(j), other->polygon()); |
| 1887 | + min_distance = std::min(min_distance, dist); |
| 1888 | + } |
| 1889 | + } |
| 1890 | + |
| 1891 | + // Check distance from each vertex of other polygon to this polygon |
| 1892 | + for (int i = 0; i < other->polygon()->num_loops(); ++i) { |
| 1893 | + const S2Loop* loop = other->polygon()->loop(i); |
| 1894 | + for (int j = 0; j < loop->num_vertices(); ++j) { |
| 1895 | + double dist = distance_point_to_polygon(loop->vertex(j), _polygon.get()); |
| 1896 | + min_distance = std::min(min_distance, dist); |
| 1897 | + } |
| 1898 | + } |
| 1899 | + |
| 1900 | + // Handle touching case: if min_distance is within tolerance, polygons are touching |
| 1901 | + if (min_distance <= TOLERANCE) { |
| 1902 | + return 0.0; |
| 1903 | + } |
| 1904 | + return min_distance; |
| 1905 | + } |
| 1906 | + case GEO_SHAPE_CIRCLE: { |
| 1907 | + return rhs->Distance(this); // Delegate to Circle's implementation |
| 1908 | + } |
| 1909 | + case GEO_SHAPE_MULTI_POLYGON: { |
| 1910 | + return rhs->Distance(this); // Delegate to MultiPolygon's implementation |
| 1911 | + } |
| 1912 | + default: |
| 1913 | + return -1.0; |
| 1914 | + } |
| 1915 | +} |
| 1916 | + |
| 1917 | +double GeoMultiPolygon::Distance(const GeoShape* rhs) const { |
| 1918 | + // rhs is guaranteed to be valid by StDistance (functions_geo.cpp) |
| 1919 | + double min_distance = std::numeric_limits<double>::max(); |
| 1920 | + |
| 1921 | + // Calculate minimum distance from any polygon to the other shape |
| 1922 | + for (const auto& polygon : _polygons) { |
| 1923 | + double dist = polygon->Distance(rhs); |
| 1924 | + if (dist >= 0) { |
| 1925 | + min_distance = std::min(min_distance, dist); |
| 1926 | + } |
| 1927 | + } |
| 1928 | + |
| 1929 | + return (min_distance == std::numeric_limits<double>::max()) ? -1.0 : min_distance; |
| 1930 | +} |
| 1931 | + |
| 1932 | +double GeoCircle::Distance(const GeoShape* rhs) const { |
| 1933 | + // Both rhs and self are guaranteed to be valid by StDistance (functions_geo.cpp) |
| 1934 | + double circle_radius = S2Earth::ToMeters(_cap->radius()); |
| 1935 | + |
| 1936 | + switch (rhs->type()) { |
| 1937 | + case GEO_SHAPE_POINT: { |
| 1938 | + return rhs->Distance(this); // Reuse Point's Distance implementation |
| 1939 | + } |
| 1940 | + case GEO_SHAPE_LINE_STRING: { |
| 1941 | + const GeoLine* line = static_cast<const GeoLine*>(rhs); |
| 1942 | + double min_distance = std::numeric_limits<double>::max(); |
| 1943 | + |
| 1944 | + // Find minimum distance from circle center to line |
| 1945 | + for (int i = 0; i < line->polyline()->num_vertices() - 1; ++i) { |
| 1946 | + double dist = distance_point_to_segment(_cap->center(), line->polyline()->vertex(i), |
| 1947 | + line->polyline()->vertex(i + 1)); |
| 1948 | + min_distance = std::min(min_distance, dist); |
| 1949 | + } |
| 1950 | + |
| 1951 | + if (min_distance <= circle_radius + TOLERANCE) { |
| 1952 | + return 0.0; |
| 1953 | + } |
| 1954 | + return min_distance - circle_radius; |
| 1955 | + } |
| 1956 | + case GEO_SHAPE_POLYGON: { |
| 1957 | + const GeoPolygon* polygon = static_cast<const GeoPolygon*>(rhs); |
| 1958 | + |
| 1959 | + // If center is inside polygon, distance is 0 |
| 1960 | + if (polygon->polygon()->Contains(_cap->center())) { |
| 1961 | + return 0.0; |
| 1962 | + } |
| 1963 | + |
| 1964 | + // Find minimum distance from circle center to polygon boundary |
| 1965 | + double min_distance = std::numeric_limits<double>::max(); |
| 1966 | + |
| 1967 | + for (int i = 0; i < polygon->polygon()->num_loops(); ++i) { |
| 1968 | + const S2Loop* loop = polygon->polygon()->loop(i); |
| 1969 | + for (int j = 0; j < loop->num_vertices(); ++j) { |
| 1970 | + double dist = |
| 1971 | + distance_point_to_segment(_cap->center(), loop->vertex(j), |
| 1972 | + loop->vertex((j + 1) % loop->num_vertices())); |
| 1973 | + min_distance = std::min(min_distance, dist); |
| 1974 | + } |
| 1975 | + } |
| 1976 | + |
| 1977 | + if (min_distance <= circle_radius + TOLERANCE) { |
| 1978 | + return 0.0; |
| 1979 | + } |
| 1980 | + return min_distance - circle_radius; |
| 1981 | + } |
| 1982 | + case GEO_SHAPE_CIRCLE: { |
| 1983 | + const GeoCircle* other = static_cast<const GeoCircle*>(rhs); |
| 1984 | + double other_radius = S2Earth::ToMeters(other->circle()->radius()); |
| 1985 | + S2LatLng this_center_ll = S2LatLng(_cap->center()); |
| 1986 | + S2LatLng other_center_ll = S2LatLng(other->circle()->center()); |
| 1987 | + double dist_centers = S2Earth::GetDistanceMeters(this_center_ll, other_center_ll); |
| 1988 | + |
| 1989 | + // Distance between circles is distance between centers minus sum of radii |
| 1990 | + double sum_radii = circle_radius + other_radius; |
| 1991 | + if (dist_centers <= sum_radii + TOLERANCE) { |
| 1992 | + return 0.0; // Circles intersect or touch |
| 1993 | + } |
| 1994 | + return dist_centers - sum_radii; |
| 1995 | + } |
| 1996 | + case GEO_SHAPE_MULTI_POLYGON: { |
| 1997 | + return rhs->Distance(this); // Delegate to MultiPolygon's implementation |
| 1998 | + } |
| 1999 | + default: |
| 2000 | + return -1.0; |
| 2001 | + } |
| 2002 | +} |
| 2003 | + |
1633 | 2004 | #include "common/compile_check_avoid_end.h" |
1634 | 2005 | } // namespace doris |
0 commit comments