/src/h3/src/h3lib/lib/bbox.c
Line | Count | Source (jump to first uncovered line) |
1 | | /* |
2 | | * Copyright 2016-2021 Uber Technologies, Inc. |
3 | | * |
4 | | * Licensed under the Apache License, Version 2.0 (the "License"); |
5 | | * you may not use this file except in compliance with the License. |
6 | | * You may obtain a copy of the License at |
7 | | * |
8 | | * http://www.apache.org/licenses/LICENSE-2.0 |
9 | | * |
10 | | * Unless required by applicable law or agreed to in writing, software |
11 | | * distributed under the License is distributed on an "AS IS" BASIS, |
12 | | * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
13 | | * See the License for the specific language governing permissions and |
14 | | * limitations under the License. |
15 | | */ |
16 | | /** @file bbox.c |
17 | | * @brief Geographic bounding box functions |
18 | | */ |
19 | | |
20 | | #include "bbox.h" |
21 | | |
22 | | #include <float.h> |
23 | | #include <math.h> |
24 | | #include <stdbool.h> |
25 | | |
26 | | #include "constants.h" |
27 | | #include "h3Index.h" |
28 | | #include "latLng.h" |
29 | | |
30 | | /** |
31 | | * Whether the given bounding box crosses the antimeridian |
32 | | * @param bbox Bounding box to inspect |
33 | | * @return is transmeridian |
34 | | */ |
35 | 54.6M | bool bboxIsTransmeridian(const BBox *bbox) { return bbox->east < bbox->west; } |
36 | | |
37 | | /** |
38 | | * Get the center of a bounding box |
39 | | * @param bbox Input bounding box |
40 | | * @param center Output center coordinate |
41 | | */ |
42 | 0 | void bboxCenter(const BBox *bbox, LatLng *center) { |
43 | 0 | center->lat = (bbox->north + bbox->south) / 2.0; |
44 | | // If the bbox crosses the antimeridian, shift east 360 degrees |
45 | 0 | double east = bboxIsTransmeridian(bbox) ? bbox->east + M_2PI : bbox->east; |
46 | 0 | center->lng = constrainLng((east + bbox->west) / 2.0); |
47 | 0 | } |
48 | | |
49 | | /** |
50 | | * Whether the bounding box contains a given point |
51 | | * @param bbox Bounding box |
52 | | * @param point Point to test |
53 | | * @return Whether the point is contained |
54 | | */ |
55 | 45.5M | bool bboxContains(const BBox *bbox, const LatLng *point) { |
56 | 45.5M | return point->lat >= bbox->south && point->lat <= bbox->north && |
57 | 45.5M | (bboxIsTransmeridian(bbox) ? |
58 | | // transmeridian case |
59 | 28.7M | (point->lng >= bbox->west || point->lng <= bbox->east) |
60 | 30.0M | : |
61 | | // standard case |
62 | 30.0M | (point->lng >= bbox->west && point->lng <= bbox->east)); |
63 | 45.5M | } |
64 | | |
65 | | /** |
66 | | * Whether two bounding boxes are strictly equal |
67 | | * @param b1 Bounding box 1 |
68 | | * @param b2 Bounding box 2 |
69 | | * @return Whether the boxes are equal |
70 | | */ |
71 | 0 | bool bboxEquals(const BBox *b1, const BBox *b2) { |
72 | 0 | return b1->north == b2->north && b1->south == b2->south && |
73 | 0 | b1->east == b2->east && b1->west == b2->west; |
74 | 0 | } |
75 | | |
76 | | /** |
77 | | * _hexRadiusKm returns the radius of a given hexagon in Km |
78 | | * |
79 | | * @param h3Index the index of the hexagon |
80 | | * @return the radius of the hexagon in Km |
81 | | */ |
82 | 200k | double _hexRadiusKm(H3Index h3Index) { |
83 | | // There is probably a cheaper way to determine the radius of a |
84 | | // hexagon, but this way is conceptually simple |
85 | 200k | LatLng h3Center; |
86 | 200k | CellBoundary h3Boundary; |
87 | 200k | H3_EXPORT(cellToLatLng)(h3Index, &h3Center); |
88 | 200k | H3_EXPORT(cellToBoundary)(h3Index, &h3Boundary); |
89 | 200k | return H3_EXPORT(greatCircleDistanceKm)(&h3Center, h3Boundary.verts); |
90 | 200k | } |
91 | | |
92 | | /** |
93 | | * bboxHexEstimate returns an estimated number of hexagons that fit |
94 | | * within the cartesian-projected bounding box |
95 | | * |
96 | | * @param bbox the bounding box to estimate the hexagon fill level |
97 | | * @param res the resolution of the H3 hexagons to fill the bounding box |
98 | | * @param out the estimated number of hexagons to fill the bounding box |
99 | | * @return E_SUCCESS (0) on success, or another value otherwise. |
100 | | */ |
101 | 1.94k | H3Error bboxHexEstimate(const BBox *bbox, int res, int64_t *out) { |
102 | | // Get the area of the pentagon as the maximally-distorted area possible |
103 | 1.94k | H3Index pentagons[12] = {0}; |
104 | 1.94k | H3Error pentagonsErr = H3_EXPORT(getPentagons)(res, pentagons); |
105 | 1.94k | if (pentagonsErr) { |
106 | 12 | return pentagonsErr; |
107 | 12 | } |
108 | 1.93k | double pentagonRadiusKm = _hexRadiusKm(pentagons[0]); |
109 | | // Area of a regular hexagon is 3/2*sqrt(3) * r * r |
110 | | // The pentagon has the most distortion (smallest edges) and shares its |
111 | | // edges with hexagons, so the most-distorted hexagons have this area, |
112 | | // shrunk by 20% off chance that the bounding box perfectly bounds a |
113 | | // pentagon. |
114 | 1.93k | double pentagonAreaKm2 = |
115 | 1.93k | 0.8 * (2.59807621135 * pentagonRadiusKm * pentagonRadiusKm); |
116 | | |
117 | | // Then get the area of the bounding box of the geoloop in question |
118 | 1.93k | LatLng p1, p2; |
119 | 1.93k | p1.lat = bbox->north; |
120 | 1.93k | p1.lng = bbox->east; |
121 | 1.93k | p2.lat = bbox->south; |
122 | 1.93k | p2.lng = bbox->west; |
123 | 1.93k | double d = H3_EXPORT(greatCircleDistanceKm)(&p1, &p2); |
124 | 1.93k | double lngDiff = fabs(p1.lng - p2.lng); |
125 | 1.93k | double latDiff = fabs(p1.lat - p2.lat); |
126 | 1.93k | if (lngDiff == 0 || latDiff == 0) { |
127 | 16 | return E_FAILED; |
128 | 16 | } |
129 | 1.92k | double length = fmax(lngDiff, latDiff); |
130 | 1.92k | double width = fmin(lngDiff, latDiff); |
131 | 1.92k | double ratio = length / width; |
132 | | // Derived constant based on: https://math.stackexchange.com/a/1921940 |
133 | | // Clamped to 3 as higher values tend to rapidly drag the estimate to zero. |
134 | 1.92k | double a = d * d / fmin(3.0, ratio); |
135 | | |
136 | | // Divide the two to get an estimate of the number of hexagons needed |
137 | 1.92k | double estimateDouble = ceil(a / pentagonAreaKm2); |
138 | 1.92k | if (!isfinite(estimateDouble)) { |
139 | 4 | return E_FAILED; |
140 | 4 | } |
141 | 1.91k | int64_t estimate = (int64_t)estimateDouble; |
142 | 1.91k | if (estimate == 0) estimate = 1; |
143 | 1.91k | *out = estimate; |
144 | 1.91k | return E_SUCCESS; |
145 | 1.92k | } |
146 | | |
147 | | /** |
148 | | * lineHexEstimate returns an estimated number of hexagons that trace |
149 | | * the cartesian-projected line |
150 | | * |
151 | | * @param origin the origin coordinates |
152 | | * @param destination the destination coordinates |
153 | | * @param res the resolution of the H3 hexagons to trace the line |
154 | | * @param out Out parameter for the estimated number of hexagons required to |
155 | | * trace the line |
156 | | * @return E_SUCCESS (0) on success or another value otherwise. |
157 | | */ |
158 | | H3Error lineHexEstimate(const LatLng *origin, const LatLng *destination, |
159 | 198k | int res, int64_t *out) { |
160 | | // Get the area of the pentagon as the maximally-distorted area possible |
161 | 198k | H3Index pentagons[12] = {0}; |
162 | 198k | H3Error pentagonsErr = H3_EXPORT(getPentagons)(res, pentagons); |
163 | 198k | if (pentagonsErr) { |
164 | 0 | return pentagonsErr; |
165 | 0 | } |
166 | 198k | double pentagonRadiusKm = _hexRadiusKm(pentagons[0]); |
167 | | |
168 | 198k | double dist = H3_EXPORT(greatCircleDistanceKm)(origin, destination); |
169 | 198k | double distCeil = ceil(dist / (2 * pentagonRadiusKm)); |
170 | 198k | if (!isfinite(distCeil)) { |
171 | 148 | return E_FAILED; |
172 | 148 | } |
173 | 198k | int64_t estimate = (int64_t)distCeil; |
174 | 198k | if (estimate == 0) estimate = 1; |
175 | 198k | *out = estimate; |
176 | 198k | return E_SUCCESS; |
177 | 198k | } |