1/*
2 * Copyright 2018 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 Geofence 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#include "bbox.h"
29#include "constants.h"
30#include "geoCoord.h"
31#include "h3api.h"
32#include "linkedGeo.h"
33#include "polygon.h"
34
35#ifndef TYPE
36#error "TYPE must be defined before including this header"
37#endif
38
39#ifndef IS_EMPTY
40#error "IS_EMPTY must be defined before including this header"
41#endif
42
43#ifndef INIT_ITERATION
44#error "INIT_ITERATION must be defined before including this header"
45#endif
46
47#ifndef ITERATE
48#error "ITERATE must be defined before including this header"
49#endif
50
51#define LOOP_ALGO_XTJOIN(a, b) a##b
52#define LOOP_ALGO_TJOIN(a, b) LOOP_ALGO_XTJOIN(a, b)
53#define GENERIC_LOOP_ALGO(func) LOOP_ALGO_TJOIN(func, TYPE)
54
55/** Macro: Normalize longitude, dealing with transmeridian arcs */
56#define NORMALIZE_LON(lon, isTransmeridian) \
57 (isTransmeridian && lon < 0 ? lon + (double)M_2PI : lon)
58
59/**
60 * pointInside is the core loop of the point-in-poly algorithm
61 * @param loop The loop to check
62 * @param bbox The bbox for the loop being tested
63 * @param coord The coordinate to check
64 * @return Whether the point is contained
65 */
66bool GENERIC_LOOP_ALGO(pointInside)(const TYPE* loop, const BBox* bbox,
67 const GeoCoord* coord) {
68 // fail fast if we're outside the bounding box
69 if (!bboxContains(bbox, coord)) {
70 return false;
71 }
72 bool isTransmeridian = bboxIsTransmeridian(bbox);
73 bool contains = false;
74
75 double lat = coord->lat;
76 double lng = NORMALIZE_LON(coord->lon, isTransmeridian);
77
78 GeoCoord a;
79 GeoCoord b;
80
81 INIT_ITERATION;
82
83 while (true) {
84 ITERATE(loop, a, b);
85
86 // Ray casting algo requires the second point to always be higher
87 // than the first, so swap if needed
88 if (a.lat > b.lat) {
89 GeoCoord tmp = a;
90 a = b;
91 b = tmp;
92 }
93
94 // If we're totally above or below the latitude ranges, the test
95 // ray cannot intersect the line segment, so let's move on
96 if (lat < a.lat || lat > b.lat) {
97 continue;
98 }
99
100 double aLng = NORMALIZE_LON(a.lon, isTransmeridian);
101 double bLng = NORMALIZE_LON(b.lon, isTransmeridian);
102
103 // Rays are cast in the longitudinal direction, in case a point
104 // exactly matches, to decide tiebreakers, bias westerly
105 if (aLng == lng || bLng == lng) {
106 lng -= DBL_EPSILON;
107 }
108
109 // For the latitude of the point, compute the longitude of the
110 // point that lies on the line segment defined by a and b
111 // This is done by computing the percent above a the lat is,
112 // and traversing the same percent in the longitudinal direction
113 // of a to b
114 double ratio = (lat - a.lat) / (b.lat - a.lat);
115 double testLng =
116 NORMALIZE_LON(aLng + (bLng - aLng) * ratio, isTransmeridian);
117
118 // Intersection of the ray
119 if (testLng > lng) {
120 contains = !contains;
121 }
122 }
123
124 return contains;
125}
126
127/**
128 * Create a bounding box from a simple polygon loop.
129 * Known limitations:
130 * - Does not support polygons with two adjacent points > 180 degrees of
131 * longitude apart. These will be interpreted as crossing the antimeridian.
132 * - Does not currently support polygons containing a pole.
133 * @param loop Loop of coordinates
134 * @param bbox Output bbox
135 */
136void GENERIC_LOOP_ALGO(bboxFrom)(const TYPE* loop, BBox* bbox) {
137 // Early exit if there are no vertices
138 if (IS_EMPTY(loop)) {
139 *bbox = (BBox){0};
140 return;
141 }
142
143 bbox->south = DBL_MAX;
144 bbox->west = DBL_MAX;
145 bbox->north = -DBL_MAX;
146 bbox->east = -DBL_MAX;
147 double minPosLon = DBL_MAX;
148 double maxNegLon = -DBL_MAX;
149 bool isTransmeridian = false;
150
151 double lat;
152 double lon;
153 GeoCoord coord;
154 GeoCoord next;
155
156 INIT_ITERATION;
157
158 while (true) {
159 ITERATE(loop, coord, next);
160
161 lat = coord.lat;
162 lon = coord.lon;
163 if (lat < bbox->south) bbox->south = lat;
164 if (lon < bbox->west) bbox->west = lon;
165 if (lat > bbox->north) bbox->north = lat;
166 if (lon > bbox->east) bbox->east = lon;
167 // Save the min positive and max negative longitude for
168 // use in the transmeridian case
169 if (lon > 0 && lon < minPosLon) minPosLon = lon;
170 if (lon < 0 && lon > maxNegLon) maxNegLon = lon;
171 // check for arcs > 180 degrees longitude, flagging as transmeridian
172 if (fabs(lon - next.lon) > M_PI) {
173 isTransmeridian = true;
174 }
175 }
176 // Swap east and west if transmeridian
177 if (isTransmeridian) {
178 bbox->east = maxNegLon;
179 bbox->west = minPosLon;
180 }
181}
182
183/**
184 * Whether the winding order of a given loop is clockwise, with normalization
185 * for loops crossing the antimeridian.
186 * @param loop The loop to check
187 * @param isTransmeridian Whether the loop crosses the antimeridian
188 * @return Whether the loop is clockwise
189 */
190static bool GENERIC_LOOP_ALGO(isClockwiseNormalized)(const TYPE* loop,
191 bool isTransmeridian) {
192 double sum = 0;
193 GeoCoord a;
194 GeoCoord b;
195
196 INIT_ITERATION;
197 while (true) {
198 ITERATE(loop, a, b);
199 // If we identify a transmeridian arc (> 180 degrees longitude),
200 // start over with the transmeridian flag set
201 if (!isTransmeridian && fabs(a.lon - b.lon) > M_PI) {
202 return GENERIC_LOOP_ALGO(isClockwiseNormalized)(loop, true);
203 }
204 sum += ((NORMALIZE_LON(b.lon, isTransmeridian) -
205 NORMALIZE_LON(a.lon, isTransmeridian)) *
206 (b.lat + a.lat));
207 }
208
209 return sum > 0;
210}
211
212/**
213 * Whether the winding order of a given loop is clockwise. In GeoJSON,
214 * clockwise loops are always inner loops (holes).
215 * @param loop The loop to check
216 * @return Whether the loop is clockwise
217 */
218bool GENERIC_LOOP_ALGO(isClockwise)(const TYPE* loop) {
219 return GENERIC_LOOP_ALGO(isClockwiseNormalized)(loop, false);
220}
221