/src/h3/src/h3lib/include/polygonAlgos.h
Line | Count | Source (jump to first uncovered line) |
1 | | /* |
2 | | * Copyright 2018, 2020-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 |
17 | | * @brief Include file for poylgon algorithms. This includes the core |
18 | | * logic for algorithms acting over loops of coordinates, |
19 | | * allowing them to be reused for both GeoLoop and |
20 | | * LinkegGeoLoop structures. This file is intended to be |
21 | | * included inline in a file that defines the type-specific |
22 | | * macros required for iteration. |
23 | | */ |
24 | | |
25 | | #include <float.h> |
26 | | #include <math.h> |
27 | | #include <stdbool.h> |
28 | | |
29 | | #include "bbox.h" |
30 | | #include "constants.h" |
31 | | #include "h3api.h" |
32 | | #include "latLng.h" |
33 | | #include "linkedGeo.h" |
34 | | #include "polygon.h" |
35 | | |
36 | | #ifndef TYPE |
37 | | #error "TYPE must be defined before including this header" |
38 | | #endif |
39 | | |
40 | | #ifndef IS_EMPTY |
41 | | #error "IS_EMPTY must be defined before including this header" |
42 | | #endif |
43 | | |
44 | | #ifndef INIT_ITERATION |
45 | | #error "INIT_ITERATION must be defined before including this header" |
46 | | #endif |
47 | | |
48 | | #ifndef ITERATE |
49 | | #error "ITERATE must be defined before including this header" |
50 | | #endif |
51 | | |
52 | 60.8k | #define LOOP_ALGO_XTJOIN(a, b) a##b |
53 | 60.8k | #define LOOP_ALGO_TJOIN(a, b) LOOP_ALGO_XTJOIN(a, b) |
54 | 60.8k | #define GENERIC_LOOP_ALGO(func) LOOP_ALGO_TJOIN(func, TYPE) |
55 | | |
56 | | /** Macro: Normalize longitude, dealing with transmeridian arcs */ |
57 | | #define NORMALIZE_LNG(lng, isTransmeridian) \ |
58 | 605k | (isTransmeridian && lng < 0 ? lng + (double)M_2PI : lng) |
59 | | |
60 | | /** |
61 | | * pointInside is the core loop of the point-in-poly algorithm |
62 | | * @param loop The loop to check |
63 | | * @param bbox The bbox for the loop being tested |
64 | | * @param coord The coordinate to check |
65 | | * @return Whether the point is contained |
66 | | */ |
67 | | bool GENERIC_LOOP_ALGO(pointInside)(const TYPE *loop, const BBox *bbox, |
68 | 2.46M | const LatLng *coord) { |
69 | | // fail fast if we're outside the bounding box |
70 | 2.46M | if (!bboxContains(bbox, coord)) { |
71 | 2.45M | return false; |
72 | 2.45M | } |
73 | 7.54k | bool isTransmeridian = bboxIsTransmeridian(bbox); |
74 | 7.54k | bool contains = false; |
75 | | |
76 | 7.54k | double lat = coord->lat; |
77 | 7.54k | double lng = NORMALIZE_LNG(coord->lng, isTransmeridian); |
78 | | |
79 | 7.54k | LatLng a; |
80 | 7.54k | LatLng b; |
81 | | |
82 | 7.54k | INIT_ITERATION; |
83 | | |
84 | 191k | while (true) { |
85 | 191k | ITERATE(loop, a, b); |
86 | | |
87 | | // Ray casting algo requires the second point to always be higher |
88 | | // than the first, so swap if needed |
89 | 183k | if (a.lat > b.lat) { |
90 | 93.7k | LatLng tmp = a; |
91 | 93.7k | a = b; |
92 | 93.7k | b = tmp; |
93 | 93.7k | } |
94 | | |
95 | | // If the latitude matches exactly, we'll hit an edge case where |
96 | | // the ray passes through the vertex twice on successive segment |
97 | | // checks. To avoid this, adjust the latiude northward if needed. |
98 | | // |
99 | | // NOTE: This currently means that a point at the north pole cannot |
100 | | // be contained in any polygon. This is acceptable in current usage, |
101 | | // because the point we test in this function at present is always |
102 | | // a cell center or vertex, and no cell has a center or vertex on the |
103 | | // north pole. If we need to expand this algo to more generic uses we |
104 | | // might need to handle this edge case. |
105 | 183k | if (lat == a.lat || lat == b.lat) { |
106 | 2.13k | lat += DBL_EPSILON; |
107 | 2.13k | } |
108 | | |
109 | | // If we're totally above or below the latitude ranges, the test |
110 | | // ray cannot intersect the line segment, so let's move on |
111 | 183k | if (lat < a.lat || lat > b.lat) { |
112 | 162k | continue; |
113 | 162k | } |
114 | | |
115 | 21.1k | double aLng = NORMALIZE_LNG(a.lng, isTransmeridian); |
116 | 21.1k | double bLng = NORMALIZE_LNG(b.lng, isTransmeridian); |
117 | | |
118 | | // Rays are cast in the longitudinal direction, in case a point |
119 | | // exactly matches, to decide tiebreakers, bias westerly |
120 | 21.1k | if (aLng == lng || bLng == lng) { |
121 | 1.83k | lng -= DBL_EPSILON; |
122 | 1.83k | } |
123 | | |
124 | | // For the latitude of the point, compute the longitude of the |
125 | | // point that lies on the line segment defined by a and b |
126 | | // This is done by computing the percent above a the lat is, |
127 | | // and traversing the same percent in the longitudinal direction |
128 | | // of a to b |
129 | 21.1k | double ratio = (lat - a.lat) / (b.lat - a.lat); |
130 | 21.1k | double testLng = |
131 | 21.1k | NORMALIZE_LNG(aLng + (bLng - aLng) * ratio, isTransmeridian); |
132 | | |
133 | | // Intersection of the ray |
134 | 21.1k | if (testLng > lng) { |
135 | 10.4k | contains = !contains; |
136 | 10.4k | } |
137 | 21.1k | } |
138 | | |
139 | 7.54k | return contains; |
140 | 2.46M | } Unexecuted instantiation: pointInsideGeoLoop Line | Count | Source | 68 | 2.46M | const LatLng *coord) { | 69 | | // fail fast if we're outside the bounding box | 70 | 2.46M | if (!bboxContains(bbox, coord)) { | 71 | 2.45M | return false; | 72 | 2.45M | } | 73 | 7.54k | bool isTransmeridian = bboxIsTransmeridian(bbox); | 74 | 7.54k | bool contains = false; | 75 | | | 76 | 7.54k | double lat = coord->lat; | 77 | 7.54k | double lng = NORMALIZE_LNG(coord->lng, isTransmeridian); | 78 | | | 79 | 7.54k | LatLng a; | 80 | 7.54k | LatLng b; | 81 | | | 82 | 7.54k | INIT_ITERATION; | 83 | | | 84 | 191k | while (true) { | 85 | 191k | ITERATE(loop, a, b); | 86 | | | 87 | | // Ray casting algo requires the second point to always be higher | 88 | | // than the first, so swap if needed | 89 | 183k | if (a.lat > b.lat) { | 90 | 93.7k | LatLng tmp = a; | 91 | 93.7k | a = b; | 92 | 93.7k | b = tmp; | 93 | 93.7k | } | 94 | | | 95 | | // If the latitude matches exactly, we'll hit an edge case where | 96 | | // the ray passes through the vertex twice on successive segment | 97 | | // checks. To avoid this, adjust the latiude northward if needed. | 98 | | // | 99 | | // NOTE: This currently means that a point at the north pole cannot | 100 | | // be contained in any polygon. This is acceptable in current usage, | 101 | | // because the point we test in this function at present is always | 102 | | // a cell center or vertex, and no cell has a center or vertex on the | 103 | | // north pole. If we need to expand this algo to more generic uses we | 104 | | // might need to handle this edge case. | 105 | 183k | if (lat == a.lat || lat == b.lat) { | 106 | 2.13k | lat += DBL_EPSILON; | 107 | 2.13k | } | 108 | | | 109 | | // If we're totally above or below the latitude ranges, the test | 110 | | // ray cannot intersect the line segment, so let's move on | 111 | 183k | if (lat < a.lat || lat > b.lat) { | 112 | 162k | continue; | 113 | 162k | } | 114 | | | 115 | 21.1k | double aLng = NORMALIZE_LNG(a.lng, isTransmeridian); | 116 | 21.1k | double bLng = NORMALIZE_LNG(b.lng, isTransmeridian); | 117 | | | 118 | | // Rays are cast in the longitudinal direction, in case a point | 119 | | // exactly matches, to decide tiebreakers, bias westerly | 120 | 21.1k | if (aLng == lng || bLng == lng) { | 121 | 1.83k | lng -= DBL_EPSILON; | 122 | 1.83k | } | 123 | | | 124 | | // For the latitude of the point, compute the longitude of the | 125 | | // point that lies on the line segment defined by a and b | 126 | | // This is done by computing the percent above a the lat is, | 127 | | // and traversing the same percent in the longitudinal direction | 128 | | // of a to b | 129 | 21.1k | double ratio = (lat - a.lat) / (b.lat - a.lat); | 130 | 21.1k | double testLng = | 131 | 21.1k | NORMALIZE_LNG(aLng + (bLng - aLng) * ratio, isTransmeridian); | 132 | | | 133 | | // Intersection of the ray | 134 | 21.1k | if (testLng > lng) { | 135 | 10.4k | contains = !contains; | 136 | 10.4k | } | 137 | 21.1k | } | 138 | | | 139 | 7.54k | return contains; | 140 | 2.46M | } |
|
141 | | |
142 | | /** |
143 | | * Create a bounding box from a simple polygon loop. |
144 | | * Known limitations: |
145 | | * - Does not support polygons with two adjacent points > 180 degrees of |
146 | | * longitude apart. These will be interpreted as crossing the antimeridian. |
147 | | * - Does not currently support polygons containing a pole. |
148 | | * @param loop Loop of coordinates |
149 | | * @param bbox Output bbox |
150 | | */ |
151 | 58.7k | void GENERIC_LOOP_ALGO(bboxFrom)(const TYPE *loop, BBox *bbox) { |
152 | | // Early exit if there are no vertices |
153 | 58.7k | if (IS_EMPTY(loop)) { |
154 | 0 | *bbox = (BBox){0}; |
155 | 0 | return; |
156 | 0 | } |
157 | | |
158 | 58.7k | bbox->south = DBL_MAX; |
159 | 58.7k | bbox->west = DBL_MAX; |
160 | 58.7k | bbox->north = -DBL_MAX; |
161 | 58.7k | bbox->east = -DBL_MAX; |
162 | 58.7k | double minPosLng = DBL_MAX; |
163 | 58.7k | double maxNegLng = -DBL_MAX; |
164 | 58.7k | bool isTransmeridian = false; |
165 | | |
166 | 58.7k | double lat; |
167 | 58.7k | double lng; |
168 | 58.7k | LatLng coord; |
169 | 58.7k | LatLng next; |
170 | | |
171 | 58.7k | INIT_ITERATION; |
172 | | |
173 | 314k | while (true) { |
174 | 314k | ITERATE(loop, coord, next); |
175 | | |
176 | 256k | lat = coord.lat; |
177 | 256k | lng = coord.lng; |
178 | 256k | if (lat < bbox->south) bbox->south = lat; |
179 | 256k | if (lng < bbox->west) bbox->west = lng; |
180 | 256k | if (lat > bbox->north) bbox->north = lat; |
181 | 256k | if (lng > bbox->east) bbox->east = lng; |
182 | | // Save the min positive and max negative longitude for |
183 | | // use in the transmeridian case |
184 | 256k | if (lng > 0 && lng < minPosLng) minPosLng = lng; |
185 | 256k | if (lng < 0 && lng > maxNegLng) maxNegLng = lng; |
186 | | // check for arcs > 180 degrees longitude, flagging as transmeridian |
187 | 256k | if (fabs(lng - next.lng) > M_PI) { |
188 | 1.00k | isTransmeridian = true; |
189 | 1.00k | } |
190 | 256k | } |
191 | | // Swap east and west if transmeridian |
192 | 58.7k | if (isTransmeridian) { |
193 | 478 | bbox->east = maxNegLng; |
194 | 478 | bbox->west = minPosLng; |
195 | 478 | } |
196 | 58.7k | } Unexecuted instantiation: bboxFromGeoLoop Line | Count | Source | 151 | 58.7k | void GENERIC_LOOP_ALGO(bboxFrom)(const TYPE *loop, BBox *bbox) { | 152 | | // Early exit if there are no vertices | 153 | 58.7k | if (IS_EMPTY(loop)) { | 154 | 0 | *bbox = (BBox){0}; | 155 | 0 | return; | 156 | 0 | } | 157 | | | 158 | 58.7k | bbox->south = DBL_MAX; | 159 | 58.7k | bbox->west = DBL_MAX; | 160 | 58.7k | bbox->north = -DBL_MAX; | 161 | 58.7k | bbox->east = -DBL_MAX; | 162 | 58.7k | double minPosLng = DBL_MAX; | 163 | 58.7k | double maxNegLng = -DBL_MAX; | 164 | 58.7k | bool isTransmeridian = false; | 165 | | | 166 | 58.7k | double lat; | 167 | 58.7k | double lng; | 168 | 58.7k | LatLng coord; | 169 | 58.7k | LatLng next; | 170 | | | 171 | 58.7k | INIT_ITERATION; | 172 | | | 173 | 314k | while (true) { | 174 | 314k | ITERATE(loop, coord, next); | 175 | | | 176 | 256k | lat = coord.lat; | 177 | 256k | lng = coord.lng; | 178 | 256k | if (lat < bbox->south) bbox->south = lat; | 179 | 256k | if (lng < bbox->west) bbox->west = lng; | 180 | 256k | if (lat > bbox->north) bbox->north = lat; | 181 | 256k | if (lng > bbox->east) bbox->east = lng; | 182 | | // Save the min positive and max negative longitude for | 183 | | // use in the transmeridian case | 184 | 256k | if (lng > 0 && lng < minPosLng) minPosLng = lng; | 185 | 256k | if (lng < 0 && lng > maxNegLng) maxNegLng = lng; | 186 | | // check for arcs > 180 degrees longitude, flagging as transmeridian | 187 | 256k | if (fabs(lng - next.lng) > M_PI) { | 188 | 1.00k | isTransmeridian = true; | 189 | 1.00k | } | 190 | 256k | } | 191 | | // Swap east and west if transmeridian | 192 | 58.7k | if (isTransmeridian) { | 193 | 478 | bbox->east = maxNegLng; | 194 | 478 | bbox->west = minPosLng; | 195 | 478 | } | 196 | 58.7k | } |
|
197 | | |
198 | | /** |
199 | | * Whether the winding order of a given loop is clockwise, with normalization |
200 | | * for loops crossing the antimeridian. |
201 | | * @param loop The loop to check |
202 | | * @param isTransmeridian Whether the loop crosses the antimeridian |
203 | | * @return Whether the loop is clockwise |
204 | | */ |
205 | | static bool GENERIC_LOOP_ALGO(isClockwiseNormalized)(const TYPE *loop, |
206 | 60.8k | bool isTransmeridian) { |
207 | 60.8k | double sum = 0; |
208 | 60.8k | LatLng a; |
209 | 60.8k | LatLng b; |
210 | | |
211 | 60.8k | INIT_ITERATION; |
212 | 328k | while (true) { |
213 | 328k | ITERATE(loop, a, b); |
214 | | // If we identify a transmeridian arc (> 180 degrees longitude), |
215 | | // start over with the transmeridian flag set |
216 | 268k | if (!isTransmeridian && fabs(a.lng - b.lng) > M_PI) { |
217 | 709 | return GENERIC_LOOP_ALGO(isClockwiseNormalized)(loop, true); |
218 | 709 | } |
219 | 267k | sum += ((NORMALIZE_LNG(b.lng, isTransmeridian) - |
220 | 267k | NORMALIZE_LNG(a.lng, isTransmeridian)) * |
221 | 267k | (b.lat + a.lat)); |
222 | 267k | } |
223 | | |
224 | 60.1k | return sum > 0; |
225 | 60.8k | } Unexecuted instantiation: polygon.c:isClockwiseNormalizedGeoLoop linkedGeo.c:isClockwiseNormalizedLinkedGeoLoop Line | Count | Source | 206 | 60.8k | bool isTransmeridian) { | 207 | 60.8k | double sum = 0; | 208 | 60.8k | LatLng a; | 209 | 60.8k | LatLng b; | 210 | | | 211 | 60.8k | INIT_ITERATION; | 212 | 328k | while (true) { | 213 | 328k | ITERATE(loop, a, b); | 214 | | // If we identify a transmeridian arc (> 180 degrees longitude), | 215 | | // start over with the transmeridian flag set | 216 | 268k | if (!isTransmeridian && fabs(a.lng - b.lng) > M_PI) { | 217 | 709 | return GENERIC_LOOP_ALGO(isClockwiseNormalized)(loop, true); | 218 | 709 | } | 219 | 267k | sum += ((NORMALIZE_LNG(b.lng, isTransmeridian) - | 220 | 267k | NORMALIZE_LNG(a.lng, isTransmeridian)) * | 221 | 267k | (b.lat + a.lat)); | 222 | 267k | } | 223 | | | 224 | 60.1k | return sum > 0; | 225 | 60.8k | } |
|
226 | | |
227 | | /** |
228 | | * Whether the winding order of a given loop is clockwise. In GeoJSON, |
229 | | * clockwise loops are always inner loops (holes). |
230 | | * @param loop The loop to check |
231 | | * @return Whether the loop is clockwise |
232 | | */ |
233 | 60.1k | bool GENERIC_LOOP_ALGO(isClockwise)(const TYPE *loop) { |
234 | 60.1k | return GENERIC_LOOP_ALGO(isClockwiseNormalized)(loop, false); |
235 | 60.1k | } Unexecuted instantiation: isClockwiseGeoLoop Line | Count | Source | 233 | 60.1k | bool GENERIC_LOOP_ALGO(isClockwise)(const TYPE *loop) { | 234 | 60.1k | return GENERIC_LOOP_ALGO(isClockwiseNormalized)(loop, false); | 235 | 60.1k | } |
|