1/*
2 * Copyright 2016-2017 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#include <float.h>
22#include <math.h>
23#include <stdbool.h>
24#include "constants.h"
25#include "geoCoord.h"
26#include "h3Index.h"
27
28/**
29 * Whether the given bounding box crosses the antimeridian
30 * @param bbox Bounding box to inspect
31 * @return is transmeridian
32 */
33bool bboxIsTransmeridian(const BBox* bbox) { return bbox->east < bbox->west; }
34
35/**
36 * Get the center of a bounding box
37 * @param bbox Input bounding box
38 * @param center Output center coordinate
39 */
40void bboxCenter(const BBox* bbox, GeoCoord* center) {
41 center->lat = (bbox->north + bbox->south) / 2.0;
42 // If the bbox crosses the antimeridian, shift east 360 degrees
43 double east = bboxIsTransmeridian(bbox) ? bbox->east + M_2PI : bbox->east;
44 center->lon = constrainLng((east + bbox->west) / 2.0);
45}
46
47/**
48 * Whether the bounding box contains a given point
49 * @param bbox Bounding box
50 * @param point Point to test
51 * @return Whether the point is contained
52 */
53bool bboxContains(const BBox* bbox, const GeoCoord* point) {
54 return point->lat >= bbox->south && point->lat <= bbox->north &&
55 (bboxIsTransmeridian(bbox) ?
56 // transmeridian case
57 (point->lon >= bbox->west || point->lon <= bbox->east)
58 :
59 // standard case
60 (point->lon >= bbox->west && point->lon <= bbox->east));
61}
62
63/**
64 * Whether two bounding boxes are strictly equal
65 * @param b1 Bounding box 1
66 * @param b2 Bounding box 2
67 * @return Whether the boxes are equal
68 */
69bool bboxEquals(const BBox* b1, const BBox* b2) {
70 return b1->north == b2->north && b1->south == b2->south &&
71 b1->east == b2->east && b1->west == b2->west;
72}
73
74/**
75 * _hexRadiusKm returns the radius of a given hexagon in Km
76 *
77 * @param h3Index the index of the hexagon
78 * @return the radius of the hexagon in Km
79 */
80double _hexRadiusKm(H3Index h3Index) {
81 // There is probably a cheaper way to determine the radius of a
82 // hexagon, but this way is conceptually simple
83 GeoCoord h3Center;
84 GeoBoundary h3Boundary;
85 H3_EXPORT(h3ToGeo)(h3Index, &h3Center);
86 H3_EXPORT(h3ToGeoBoundary)(h3Index, &h3Boundary);
87 return _geoDistKm(&h3Center, h3Boundary.verts);
88}
89
90/**
91 * Get the radius of the bbox in hexagons - i.e. the radius of a k-ring centered
92 * on the bbox center and covering the entire bbox.
93 * @param bbox Bounding box to measure
94 * @param res Resolution of hexagons to use in measurement
95 * @return Radius in hexagons
96 */
97int bboxHexRadius(const BBox* bbox, int res) {
98 // Determine the center of the bounding box
99 GeoCoord center;
100 bboxCenter(bbox, &center);
101
102 // Use a vertex on the side closest to the equator, to ensure the longest
103 // radius in cases with significant distortion. East/west is arbitrary.
104 double lat =
105 fabs(bbox->north) > fabs(bbox->south) ? bbox->south : bbox->north;
106 GeoCoord vertex = {lat, bbox->east};
107
108 // Determine the length of the bounding box "radius" to then use
109 // as a circle on the earth that the k-rings must be greater than
110 double bboxRadiusKm = _geoDistKm(&center, &vertex);
111
112 // Determine the radius of the center hexagon
113 double centerHexRadiusKm = _hexRadiusKm(H3_EXPORT(geoToH3)(&center, res));
114
115 // The closest point along a hexagon drawn through the center points
116 // of a k-ring aggregation is exactly 1.5 radii of the hexagon. For
117 // any orientation of the GeoJSON encased in a circle defined by the
118 // bounding box radius and center, it is guaranteed to fit in this k-ring
119 // Rounded *up* to guarantee containment
120 return (int)ceil(bboxRadiusKm / (1.5 * centerHexRadiusKm));
121}
122