from __future__ import print_function, unicode_literals, division
import bisect
import math
import decimal
import heapq
class Angle(object):
def __init__(self, radians=0):
if not isinstance(radians, (long, float, int)):
raise ValueError()
self.__radians = radians
def __eq__(self, other):
return isinstance(other, self.__class__) \
and self.__radians == other.__radians
def __add__(self, other):
return self.__class__.from_radians(self.__radians + other.__radians)
def __repr__(self):
return '{}: {}'.format(self.__class__.__name__, self.__radians)
@classmethod
def from_degrees(cls, degrees):
return cls(math.radians(degrees))
@classmethod
def from_radians(cls, radians):
return cls(radians)
@property
def radians(self):
return self.__radians
@property
def degrees(self):
return math.degrees(self.__radians)
class Point(object):
def __init__(self, x, y, z):
self.__point = (x, y, z)
def __getitem__(self, index):
return self.__point[index]
def __neg__(self):
return self.__class__(-self[0], -self[1], -self[2])
def __eq__(self, other):
return isinstance(other, self.__class__) \
and self.__point == other.__point
def __ne_(self, other):
return not self.__eq__(other)
def __hash__(self):
return hash(self.__point)
def __repr__(self):
return 'Point: {}'.format(self.__point)
def __add__(self, other):
return self.__class__(self[0] + other[0],
self[1] + other[1],
self[2] + other[2])
def __sub__(self, other):
return self.__class__(self[0] - other[0],
self[1] - other[1],
self[2] - other[2])
def __mul__(self, other):
return self.__class__(self[0] * other,
self[1] * other,
self[2] * other)
def __rmul__(self, other):
return self.__class__(self[0] * other,
self[1] * other,
self[2] * other)
def abs(self):
return self.__class__(math.fabs(self.__point[0]),
math.fabs(self.__point[1]),
math.fabs(self.__point[2]))
def largest_abs_component(self):
temp = self.abs()
if temp[0] > temp[1]:
if temp[0] > temp[2]:
return 0
else:
return 2
else:
if temp[1] > temp[2]:
return 1
else:
return 2
def angle(self, other):
return math.atan2(self.cross_prod(other).norm(), self.dot_prod(other))
def cross_prod(self, other):
x, y, z = self.__point
ox, oy, oz = other.__point
return self.__class__(y * oz - z * oy,
z * ox - x * oz,
x * oy - y * ox)
def dot_prod(self, other):
x, y, z = self.__point
ox, oy, oz = other.__point
return x * ox + y * oy + z * oz
def norm2(self):
x, y, z = self.__point
return x * x + y * y + z * z
def norm(self):
return math.sqrt(self.norm2())
def normalize(self):
n = self.norm()
if n != 0:
n = 1.0 / n
return self.__class__(self[0] * n, self[1] * n, self[2] * n)
class LatLon(object):
@classmethod
def from_degrees(cls, lat, lon):
return cls(math.radians(lat), math.radians(lon))
@classmethod
def from_radians(cls, lat, lon):
return cls(lat, lon)
@classmethod
def from_point(cls, point):
return cls(LatLon.latitude(point).radians,
LatLon.longitude(point).radians)
@classmethod
def from_angles(cls, lat, lon):
return cls(lat.radians, lon.radians)
@classmethod
def default(cls):
return cls(0, 0)
@classmethod
def invalid(cls):
return cls(math.pi, 2 * math.pi)
def __init__(self, lat, lon):
self.__coords = (lat, lon)
def __eq__(self, other):
return isinstance(other, LatLon) and self.__coords == other.__coords
def __hash__(self):
return hash(self.__coords)
def __repr__(self):
return 'LatLon: {},{}'.format(math.degrees(self.__coords[0]),
math.degrees(self.__coords[1]))
def __add__(self, other):
return self.__class__(self.lat().radians + other.lat().radians,
self.lon().radians + other.lon().radians)
def __sub__(self, other):
return self.__class__(self.lat().radians - other.lat().radians,
self.lon().radians - other.lon().radians)
# only implemented so we can do scalar * LatLon
def __rmul__(self, other):
return self.__class__(other * self.lat().radians,
other * self.lon().radians)
@staticmethod
def latitude(point):
return Angle.from_radians(math.atan2(point[2],
math.sqrt(point[0] * point[0] + point[1] * point[1])))
@staticmethod
def longitude(point):
return Angle.from_radians(math.atan2(point[1], point[0]))
def lat(self):
return Angle.from_radians(self.__coords[0])
def lon(self):
return Angle.from_radians(self.__coords[1])
def is_valid(self):
return abs(self.lat().radians) <= (math.pi / 2) \
and abs(self.lon().radians) <= math.pi
def to_point(self):
phi = self.lat().radians
theta = self.lon().radians
cosphi = math.cos(phi)
return Point(math.cos(theta) * cosphi,
math.sin(theta) * cosphi,
math.sin(phi))
def normalized(self):
return self.__class__(
max(-math.pi / 2.0, min(math.pi / 2.0, self.lat().radians)),
drem(self.lon().radians, 2 * math.pi))
def approx_equals(self, other, max_error=1e-15):
return math.fabs(self.lat().radians - other.lat().radians) \
< max_error and \
math.fabs(self.lon().radians - other.lon().radians) < max_error
def get_distance(self, other):
assert self.is_valid()
assert other.is_valid()
from_lat = self.lat().radians
to_lat = other.lat().radians
from_lon = self.lon().radians
to_lon = other.lon().radians
dlat = math.sin(0.5 * (to_lat - from_lat))
dlon = math.sin(0.5 * (to_lon - from_lon))
x = dlat * dlat + dlon * dlon * math.cos(from_lat) * math.cos(to_lat)
return Angle.from_radians(2 * math.atan2(
math.sqrt(x), math.sqrt(max(0.0, 1.0 - x))))
class Cap(object):
ROUND_UP = 1.0 + 1.0 / (1 << 52)
def __init__(self, axis=Point(1, 0, 0), height=-1):
self.__axis = axis
self.__height = height
def __repr__(self):
return '{}: {} {}'.format(self.__class__.__name__, self.__axis,
self.__height)
@classmethod
def from_axis_height(cls, axis, height):
assert is_unit_length(axis)
return cls(axis, height)
@classmethod
def from_axis_angle(cls, axis, angle):
assert is_unit_length(axis)
assert angle.radians >= 0
return cls(axis, cls.get_height_for_angle(angle.radians))
@classmethod
def get_height_for_angle(cls, radians):
assert radians >= 0
if radians >= math.pi:
return 2
d = math.sin(0.5 * radians)
return 2 * d * d
@classmethod
def from_axis_area(cls, axis, area):
assert is_unit_length(axis)
return cls(axis, area / (2 * math.pi))
@classmethod
def empty(cls):
return cls()
@classmethod
def full(cls):
return cls(Point(1, 0, 0), 2)
def height(self):
return self.__height
def axis(self):
return self.__axis
def area(self):
return 2 * math.pi * max(0.0, self.height())
def angle(self):
if self.is_empty():
return Angle.from_radians(-1)
return Angle.from_radians(2 * math.asin(math.sqrt(0.5 * self.height())))
def is_valid(self):
return is_unit_length(self.axis()) and self.height() <= 2
def is_empty(self):
return self.height() < 0
def is_full(self):
return self.height() >= 2
def get_cap_bound(self):
return self
def add_point(self, point):
assert is_unit_length(point)
if self.is_empty():
self.__axis = point
self.__height = 0
else:
dist2 = (self.axis() - point).norm2()
self.__height = max(self.__height,
self.__class__.ROUND_UP * 0.5 * dist2)
def complement(self):
if self.is_full():
height = -1
else:
height = 2 - max(self.height(), 0.0)
return self.__class__.from_axis_height(-self.axis(), height)
def contains(self, other):
if isinstance(other, self.__class__):
if self.is_full() or other.is_empty():
return True
return self.angle().radians >= self.axis().angle(other.axis()) \
+ other.angle().radians
elif isinstance(other, Point):
assert is_unit_length(other)
return (self.axis() - other).norm2() <= 2 * self.height()
elif isinstance(other, Cell):
vertices = []
for k in xrange(4):
vertices.append(other.get_vertex(k))
if not self.contains(vertices[k]):
return False
return not self.complement().intersects(other, vertices)
else:
raise NotImplementedError()
def interior_contains(self, other):
if isinstance(other, Point):
assert is_unit_length(other)
return self.is_full() \
or (self.axis() - other).norm2() < 2 * self.height()
else:
raise NotImplementedError()
def intersects(self, *args):
if len(args) == 1 and isinstance(args[0], self.__class__):
other = args[0]
if self.is_empty() or other.is_empty():
return False
return self.angle().radians + other.angle().radians \
>= self.axis().angle(other.axis())
elif len(args) == 2 and isinstance(args[0], Cell) \
and isinstance(args[1], (list, tuple)):
cell, vertices = args
if self.height() >= 1:
return False
if self.is_empty():
return False
if cell.contains(self.axis()):
return True
sin2_angle = self.height() * (2 - self.height())
for k in xrange(4):
edge = cell.get_edge_raw(k)
dot = self.axis().dot_prod(edge)
if dot > 0:
continue
if dot * dot > sin2_angle * edge.norm2():
return False
dir = edge.cross_prod(self.axis())
if dir.dot_prod(vertices[k]) < 0 \
and dir.dot_prod(vertices[(k + 1) & 3]) > 0:
return True
return False
else:
raise NotImplementedError()
def may_intersect(self, other):
vertices = []
for k in xrange(4):
vertices.append(other.get_vertex(k))
if self.contains(vertices[k]):
return True
return self.intersects(other, vertices)
def interior_intersects(self, other):
if self.height() <= 0 or other.is_empty():
return False
return self.angle().radians + other.angle().radians > \
self.axis().angle(other.axis())
def get_rect_bound(self):
if self.is_empty():
return LatLonRect.empty()
axis_ll = LatLon.from_point(self.axis())
cap_angle = self.angle().radians
all_longitudes = False
lat, lon = [], []
lon.append(-math.pi)
lon.append(math.pi)
lat.append(axis_ll.lat().radians - cap_angle)
if lat[0] <= -math.pi / 2.0:
lat[0] = -math.pi / 2.0
all_longitudes = True
lat.append(axis_ll.lat().radians + cap_angle)
if lat[1] >= math.pi / 2.0:
lat[1] = math.pi / 2.0
all_longitudes = True
if not all_longitudes:
sin_a = math.sqrt(self.height() * (2 - self.height()))
sin_c = math.cos(axis_ll.lat().radians)
if sin_a <= sin_c:
angle_a = math.asin(sin_a / sin_c)
lon[0] = drem(axis_ll.lon().radians - angle_a, 2 * math.pi)
lon[1] = drem(axis_ll.lon().radians + angle_a, 2 * math.pi)
return LatLonRect(LineInterval(*lat), SphereInterval(*lon))
def approx_equals(self, other, max_error=1e-14):
return (self.axis().angle(other.axis()) <= max_error \
and math.fabs(self.height() - other.height()) <= max_error) \
or (self.is_empty() and other.height() <= max_error) \
or (other.is_empty() and self.height() <= max_error) \
or (self.is_full() and other.height() >= 2 - max_error) \
or (other.is_full() and self.height() >= 2 - max_error)
def expanded(self, distance):
assert distance.radians >= 0
if self.is_empty():
return self.__class__.empty()
return self.__class__.from_axis_angle(self.axis(),
self.angle() + distance)
class LatLonRect(object):
def __init__(self, *args):
if len(args) == 0:
self.__lat = LineInterval.empty()
self.__lon = SphereInterval.empty()
elif isinstance(args[0], LatLon) and isinstance(args[1], LatLon):
lo, hi = args
self.__lat = LineInterval(lo.lat().radians, hi.lat().radians)
self.__lon = SphereInterval(lo.lon().radians, hi.lon().radians)
elif isinstance(args[0], LineInterval) \
and isinstance(args[1], SphereInterval):
self.__lat, self.__lon = args
else:
raise NotImplementedError()
def __eq__(self, other):
return isinstance(other, self.__class__) \
and self.lat() == other.lat() and self.lon() == other.lon()
def __ne__(self, other):
return not self.__eq__(other)
def __repr__(self):
return '{}: {}, {}'.format(self.__class__.__name__,
self.lat(), self.lon())
def lat(self):
return self.__lat
def lon(self):
return self.__lon
def lat_lo(self):
return Angle.from_radians(self.lat().lo())
def lat_hi(self):
return Angle.from_radians(self.lat().hi())
def lon_lo(self):
return Angle.from_radians(self.lon().lo())
def lon_hi(self):
return Angle.from_radians(self.lon().hi())
def lo(self):
return LatLon.from_angles(self.lat_lo(), self.lon_lo())
def hi(self):
return LatLon.from_angles(self.lat_hi(), self.lon_hi())
# Construct a rectangle of the given size centered around the given point.
# "center" needs to be normalized, but "size" does not. The latitude
# interval of the result is clamped to [-90,90] degrees, and the longitude
# interval of the result is Full() if and only if the longitude size is
# 360 degrees or more. Examples of clamping (in degrees):
#
# center=(80,170), size=(40,60) -> lat=[60,90], lng=[140,-160]
# center=(10,40), size=(210,400) -> lat=[-90,90], lng=[-180,180]
# center=(-90,180), size=(20,50) -> lat=[-90,-80], lng=[155,-155]
@classmethod
def from_center_size(cls, center, size):
return cls.from_point(center).expanded(0.5 * size)
@classmethod
def from_point(cls, p):
assert p.is_valid()
return cls(p, p)
@classmethod
def from_point_pair(cls, a, b):
assert a.is_valid()
assert b.is_valid()
return LatLonRect(LineInterval.from_point_pair(a.lat().radians,
b.lat().radians),
SphereInterval.from_point_pair(a.lon().radians,
b.lon().radians))
@classmethod
def full_lat(cls):
return LineInterval(-math.pi / 2.0, math.pi / 2.0)
@classmethod
def full_lon(cls):
return SphereInterval.full()
@classmethod
def full(cls):
return cls(cls.full_lat(), cls.full_lon())
def is_full(self):
return self.lat() == self.__class__.full_lat() and self.lon().is_full()
def is_valid(self):
return math.fabs(self.lat().lo()) <= math.pi / 2.0 \
and math.fabs(self.lat().hi()) <= math.pi / 2.0 \
and self.lon().is_valid() \
and self.lat().is_empty() == self.lon().is_empty()
@classmethod
def empty(cls):
return cls()
def get_center(self):
return LatLon.from_radians(self.lat().get_center(),
self.lon().get_center())
def get_size(self):
return LatLon.from_radians(self.lat().get_length(),
self.lon().get_length())
def get_vertex(self, k):
# Twiddle bits to return the points in CCW order (SW, SE, NE, NW)
return LatLon.from_radians(self.lat().bound(k >> 1),
self.lon().bound((k >> 1) ^ (k & 1)))
def is_empty(self):
return self.lat().is_empty()
def is_point(self):
return self.lat().lo() == self.lat().hi() \
and self.lon().lo() == self.lon().hi()
def convolve_with_cap(self, angle):
cap = Cap.from_axis_angle(Point(1, 0, 0), angle)
r = self
for k in xrange(4):
vertex_cap = Cap.from_axis_height(self.get_vertex(k).to_point(),
cap.height())
r = r.union(vertex_cap.get_rect_bound())
return r
def contains(self, *args):
if isinstance(args[0], Point):
point = args[0]
return self.contains(LatLon.from_point(point))
elif isinstance(args[0], LatLon):
ll = args[0]
assert ll.is_valid()
return self.lat().contains(ll.lat().radians) \
and self.lon().contains(ll.lon().radians)
elif isinstance(args[0], self.__class__):
other = args[0]
return self.lat().contains(other.lat()) \
and self.lon().contains(other.lon())
elif isinstance(args[0], Cell):
cell = args[0]
return self.contains(cell.get_rect_bound())
else:
raise NotImplementedError()
def interior_contains(self, *args):
if isinstance(args[0], Point):
self.interior_contains(LatLon(args[0]))
elif isinstance(args[0], LatLon):
ll = args[0]
assert ll.is_valid()
return self.lat().interior_contains(ll.lat().radians) \
and self.lon().interior_contains(ll.lon().radians)
elif isinstance(args[0], self.__class__):
other = args[0]
return self.lat().interior_contains(other.lat()) \
and self.lon().interior_contains(other.lon())
else:
raise NotImplementedError()
def may_intersect(self, cell):
return self.intersects(cell.get_rect_bound())
def intersects(self, *args):
if isinstance(args[0], self.__class__):
other = args[0]
return self.lat().intersects(other.lat()) \
and self.lon().intersects(other.lon())
elif isinstance(args[0], Cell):
cell = args[0]
if self.is_empty():
return False
if self.contains(cell.get_center_raw()):
return True
if cell.contains(self.get_center().to_point()):
return True
if not self.intersects(cell.get_rect_bound()):
return False
cell_v = []
cell_ll = []
for i in xrange(4):
cell_v.append(cell.get_vertex(i))
cell_ll.append(LatLon.from_point(cell_v[i]))
if self.contains(cell_ll[i]):
return True
if cell.contains(self.get_vertex(i).to_point()):
return True
for i in xrange(4):
edge_lon = SphereInterval.from_point_pair(
cell_ll[i].lon().radians,
cell_ll[(i + 1) & 3].lon().radians)
if not self.lon().intersects(edge_lon):
continue
a = cell_v[i]
b = cell_v[(i + 1) & 3]
if edge_lon.contains(self.lon().lo()):
if self.__class__.intersects_lon_edge(a, b,
self.lat(), self.lon().lo()):
return True
if edge_lon.contains(self.lon().hi()):
if self.__class__.intersects_lon_edge(a, b,
self.lat(), self.lon().hi()):
return True
if self.__class__.intersects_lat_edge(a, b,
self.lat().lo(), self.lon()):
return True
if self.__class__.intersects_lat_edge(a, b,
self.lat().hi(), self.lon()):
return True
return False
else:
raise NotImplementedError()
@classmethod
def intersects_lon_edge(cls, a, b, lat, lon):
return simple_crossing(a, b,
LatLon.from_radians(lat.lo(), lon).to_point(),
LatLon.from_radians(lat.hi(), lon).to_point())
@classmethod
def intersects_lat_edge(cls, a, b, lat, lon):
assert is_unit_length(a)
assert is_unit_length(b)
z = robust_cross_prod(a, b).normalize()
if z[2] < 0:
z = -z
y = robust_cross_prod(z, Point(0, 0, 1)).normalize()
x = y.cross_prod(z)
assert is_unit_length(x)
assert x[2] >= 0
sin_lat = math.sin(lat)
if math.fabs(sin_lat) >= x[2]:
return False
assert x[2] > 0
cos_theta = sin_lat / x[2]
sin_theta = math.sqrt(1 - cos_theta * cos_theta)
theta = math.atan2(sin_theta, cos_theta)
ab_theta = SphereInterval.from_point_pair(
math.atan2(a.dot_prod(y), a.dot_prod(x)),
math.atan2(b.dot_prod(y), b.dot_prod(x)))
if ab_theta.contains(theta):
isect = x * cos_theta + y * sin_theta
if lon.contains(math.atan2(isect[1], isect[0])):
return True
if ab_theta.contains(-theta):
isect = x * cos_theta - y * sin_theta
if lon.contains(math.atan2(isect[1], isect[0])):
return True
return False
def interior_intersects(self, *args):
if isinstance(args[0], self.__class__):
other = args[0]
return self.lat().interior_intersects(other.lat()) \
and self.lon().interior_intersects(other.lon())
else:
raise NotImplementedError()
def union(self, other):
return self.__class__(self.lat().union(other.lat()),
self.lon().union(other.lon()))
def intersection(self, other):
lat = self.lat().intersection(other.lat())
lon = self.lon().intersection(other.lon())
if lat.is_empty() or lon.is_empty():
return self.__class__.empty()
return self.__class__(lat, lon)
# Return a rectangle that contains all points whose latitude distance from
# this rectangle is at most margin.lat(), and whose longitude distance
# from this rectangle is at most margin.lng(). In particular, latitudes
# are clamped while longitudes are wrapped. Note that any expansion of an
# empty interval remains empty, and both components of the given margin
# must be non-negative. "margin" does not need to be normalized.
#
# NOTE: If you are trying to grow a rectangle by a certain *distance* on
# the sphere (e.g. 5km), use the ConvolveWithCap() method instead.
def expanded(self, margin):
assert margin.lat().radians > 0
assert margin.lon().radians > 0
return self.__class__(
self.lat().expanded(margin.lat().radians) \
.intersection(self.full_lat()),
self.lon().expanded(margin.lon().radians))
def approx_equals(self, other, max_error=1e-15):
return self.lat().approx_equals(other.lat(), max_error) \
and self.lon().approx_equals(other.lon(), max_error)
def get_cap_bound(self):
if self.is_empty():
return Cap.empty()
if self.lat().lo() + self.lat().hi() < 0:
pole_z = -1
pole_angle = math.pi / 2.0 + self.lat().hi()
else:
pole_z = 1
pole_angle = math.pi / 2.0 - self.lat().lo()
pole_cap = Cap.from_axis_angle(Point(0, 0, pole_z),
Angle.from_radians(pole_angle))
lon_span = self.lon().hi() - self.lon().lo()
if drem(lon_span, 2 * math.pi) >= 0:
if lon_span < 2 * math.pi:
mid_cap = Cap.from_axis_angle(self.get_center().to_point(),
Angle.from_radians(0))
for k in xrange(4):
mid_cap.add_point(self.get_vertex(k).to_point())
if mid_cap.height() < pole_cap.height():
return mid_cap
return pole_cap
# Constants for CellId
LOOKUP_BITS = 4
SWAP_MASK = 0x01
INVERT_MASK = 0x02
POS_TO_IJ = ((0, 1, 3, 2),
(0, 2, 3, 1),
(3, 2, 0, 1),
(3, 1, 0, 2))
POS_TO_ORIENTATION = [SWAP_MASK, 0, 0, INVERT_MASK | SWAP_MASK]
LOOKUP_POS = [None] * (1 << (2 * LOOKUP_BITS + 2))
LOOKUP_IJ = [None] * (1 << (2 * LOOKUP_BITS + 2))
def _init_lookup_cell(level, i, j, orig_orientation, pos, orientation):
if level == LOOKUP_BITS:
ij = (i << LOOKUP_BITS) + j
LOOKUP_POS[(ij << 2) + orig_orientation] = (pos << 2) + orientation
LOOKUP_IJ[(pos << 2) + orig_orientation] = (ij << 2) + orientation
else:
level = level + 1
i <<= 1
j <<= 1
pos <<= 2
r = POS_TO_IJ[orientation]
for index in xrange(4):
_init_lookup_cell(level, i + (r[index] >> 1),
j + (r[index] & 1), orig_orientation,
pos + index, orientation ^ POS_TO_ORIENTATION[index])
_init_lookup_cell(0, 0, 0, 0, 0, 0)
_init_lookup_cell(0, 0, 0, SWAP_MASK, 0, SWAP_MASK)
_init_lookup_cell(0, 0, 0, INVERT_MASK, 0, INVERT_MASK)
_init_lookup_cell(0, 0, 0, SWAP_MASK | INVERT_MASK, 0, SWAP_MASK | INVERT_MASK)
class CellId(object):
# projection types
LINEAR_PROJECTION = 0
TAN_PROJECTION = 1
QUADRATIC_PROJECTION = 2
# current projection used
PROJECTION = QUADRATIC_PROJECTION
FACE_BITS = 3
NUM_FACES = 6
MAX_LEVEL = 30
POS_BITS = 2 * MAX_LEVEL + 1
MAX_SIZE = 1 << MAX_LEVEL
WRAP_OFFSET = NUM_FACES << POS_BITS
def __init__(self, *args, **kwargs):
if len(args) == 0:
self.__id = long(0)
elif len(args) == 1:
# modulus to ensure wrap around
self.__id = long(args[0]) % 0xffffffffffffffff
else:
raise ValueError()
def __repr__(self):
return 'CellId: {}'.format(self.id())
def __hash__(self):
return hash(self.id())
def __cmp__(self, other):
if not isinstance(other, self.__class__):
raise NotImplementedError('Only supports comparison with same type')
else:
return self.id().__cmp__(other.id())
@classmethod
def from_lat_lon(cls, ll):
return cls.from_point(ll.to_point())
@classmethod
def from_point(cls, p):
face, u, v = xyz_to_face_uv(p)
i = cls.st_to_ij(cls.uv_to_st(u))
j = cls.st_to_ij(cls.uv_to_st(v))
return cls.from_face_ij(face, i, j)
@classmethod
def from_face_pos_level(cls, face, pos, level):
return cls((face << cls.POS_BITS) + (pos | 1)).parent(level)
@classmethod
def from_face_ij(cls, face, i, j):
n = face << (cls.POS_BITS - 1)
bits = face & SWAP_MASK
for k in xrange(7, -1, -1):
mask = (1 << LOOKUP_BITS) - 1
bits += (((i >> (k * LOOKUP_BITS)) & mask) << (LOOKUP_BITS + 2))
bits += (((j >> (k * LOOKUP_BITS)) & mask) << 2)
bits = LOOKUP_POS[bits]
n |= (bits >> 2) << (k * 2 * LOOKUP_BITS)
bits &= (SWAP_MASK | INVERT_MASK)
return cls(n * 2 + 1)
@classmethod
def from_face_ij_wrap(cls, face, i, j):
# Convert i and j to the coordinates of a leaf cell just beyond the
# boundary of this face. This prevents 32-bit overflow in the case
# of finding the neighbors of a face cell
i = max(-1, min(cls.MAX_SIZE, i))
j = max(-1, min(cls.MAX_SIZE, j))
# Find the (u,v) coordinates corresponding to the center of cell (i,j).
# For our purposes it's sufficient to always use the linear projection
# from (s,t) to (u,v): u=2*s-1 and v=2*t-1.
scale = 1.0 / cls.MAX_SIZE
u = scale * (( i << 1) + 1 - cls.MAX_SIZE)
v = scale * (( j << 1) + 1 - cls.MAX_SIZE)
# Find the leaf cell coordinates on the adjacent face, and convert
# them to a cell id at the appropriate level. We convert from (u,v)
# back to (s,t) using s=0.5*(u+1), t=0.5*(v+1).
face, u, v = xyz_to_face_uv(face_uv_to_xyz(face, u, v))
return cls.from_face_ij(face, cls.st_to_ij(0.5 * (u + 1)),
cls.st_to_ij(0.5 * (v + 1)))
@classmethod
def from_face_ij_same(cls, face, i, j, same_face):
if same_face:
return cls.from_face_ij(face, i, j)
else:
return cls.from_face_ij_wrap(face, i, j)
@classmethod
def st_to_ij(cls, s):
return max(0, min(cls.MAX_SIZE - 1, int(math.floor(cls.MAX_SIZE * s))))
@classmethod
def lsb_for_level(cls, level):
return 1 << (2 * (cls.MAX_LEVEL - level))
def parent(self, *args):
assert self.is_valid()
if len(args) == 0:
assert not self.is_face()
new_lsb = self.lsb() << 2
return self.__class__((self.id() & -new_lsb) | new_lsb)
elif len(args) == 1:
level = args[0]
assert level >= 0
assert level <= self.level()
new_lsb = self.__class__.lsb_for_level(level)
return self.__class__((self.id() & -new_lsb) | new_lsb)
def child(self, pos):
assert self.is_valid()
assert not self.is_leaf()
new_lsb = self.lsb() >> 2
return self.__class__(self.id() + (2 * pos + 1 - 4) * new_lsb)
def contains(self, other):
assert self.is_valid()
assert other.is_valid()
return other >= self.range_min() and other <= self.range_max()
def intersects(self, other):
assert self.is_valid()
assert other.is_valid()
return other.range_min() <= self.range_max() \
and other.range_max() >= self.range_min()
def is_face(self):
return (self.id() & (self.lsb_for_level(0) - 1)) == 0
def id(self):
return self.__id
def is_valid(self):
return (self.face() < self.__class__.NUM_FACES) \
and (self.lsb() & 0x1555555555555555) != 0
def lsb(self):
return self.id() & -self.id()
def face(self):
return self.id() >> self.__class__.POS_BITS
def pos(self):
# is this correct?
return self.id() & (0xffffffffffffffff >> self.__class__.FACE_BITS)
def is_leaf(self):
return (self.id() & 1) != 0
def level(self):
if self.is_leaf():
return self.__class__.MAX_LEVEL
x = (self.id() & 0xffffffff)
level = -1
if x != 0:
level += 16
else:
# is this correct?
x = ((self.id() >> 32) & 0xffffffff)
x &= -x
if x & 0x00005555:
level += 8
if x & 0x00550055:
level += 4
if x & 0x05050505:
level += 2
if x & 0x11111111:
level += 1
assert level >= 0
assert level <= self.__class__.MAX_LEVEL
return level
def child_begin(self, *args):
assert self.is_valid()
if len(args) == 0:
assert not self.is_leaf()
old_lsb = self.lsb()
return self.__class__(self.id() - old_lsb + (old_lsb >> 2))
elif len(args) == 1:
level = args[0]
assert level >= self.level()
assert level <= self.__class__.MAX_LEVEL
return self.__class__(self.id() - self.lsb() \
+ self.__class__.lsb_for_level(level))
else:
raise ValueError('no args or level arg')
def child_end(self, *args):
assert self.is_valid()
if len(args) == 0:
assert not self.is_leaf()
old_lsb = self.lsb()
return self.__class__(self.id() + old_lsb + (old_lsb >> 2))
elif len(args) == 1:
level = args[0]
assert level >= self.level()
assert level <= self.__class__.MAX_LEVEL
return self.__class__(self.id() + self.lsb() \
+ self.__class__.lsb_for_level(level))
else:
raise ValueError('no args or level arg')
def prev(self):
return self.__class__(self.id() - (self.lsb() << 1))
def next(self):
return self.__class__(self.id() + (self.lsb() << 1))
def children(self, *args):
cell_id = self.child_begin(*args)
end = self.child_end(*args)
while cell_id != end:
yield cell_id
cell_id = cell_id.next()
def range_min(self):
return self.__class__(self.id() - (self.lsb() - 1))
def range_max(self):
return self.__class__(self.id() + (self.lsb() - 1))
@classmethod
def begin(cls, level):
return cls.from_face_pos_level(0, 0, 0).child_begin(level)
@classmethod
def end(cls, level):
return cls.from_face_pos_level(5, 0, 0).child_end(level)
@classmethod
def none(cls):
return cls()
def prev_wrap(self):
assert self.is_valid()
p = self.prev()
if p.id() < self.__class__.WRAP_OFFSET:
return p
else:
return self.__class__(p.id() + self.__class__.WRAP_OFFSET)
def next_wrap(self):
assert self.is_valid()
n = self.next()
if n.id() < self.__class__.WRAP_OFFSET:
return n
else:
return self.__class__(n.id() - self.__class__.WRAP_OFFSET)
def advance_wrap(self, steps):
assert self.is_valid()
if steps == 0:
return self
step_shift = 2 * (self.__class__.MAX_LEVEL - self.level()) + 1
if steps < 0:
min_steps = -(self.id() >> step_shift)
if steps < min_steps:
step_wrap = self.__class__.WRAP_OFFSET >> step_shift
# cannot use steps %= step_wrap as Python % different to C++
steps = long(math.fmod(steps, step_wrap))
if steps < min_steps:
steps += step_wrap
else:
max_steps = (self.__class__.WRAP_OFFSET - self.id()) >> step_shift
if steps > max_steps:
step_wrap = self.__class__.WRAP_OFFSET >> step_shift
# cannot use steps %= step_wrap as Python % different to C++
steps = long(math.fmod(steps, step_wrap))
if steps > max_steps:
steps -= step_wrap
return self.__class__(self.id() + (steps << step_shift))
def advance(self, steps):
if steps == 0:
return self
step_shift = 2 * (self.__class__.MAX_LEVEL - self.level()) + 1
if steps < 0:
min_steps = -(self.id() >> step_shift)
if steps < min_steps:
steps = min_steps
else:
max_steps = (self.__class__.WRAP_OFFSET + self.lsb() - self.id()) \
>> step_shift
if steps > max_steps:
steps = max_steps
return self.__class__(self.id() + (steps << step_shift))
# Return the S2LatLng corresponding to the center of the given cell.
def to_lat_lon(self):
return LatLon.from_point(self.to_point_raw())
def to_point_raw(self):
face, si, ti = self.get_center_si_ti()
return face_uv_to_xyz(face,
self.__class__.st_to_uv((0.5 / self.__class__.MAX_SIZE) * si),
self.__class__.st_to_uv((0.5 / self.__class__.MAX_SIZE) * ti))
def to_point(self):
return self.to_point_raw().normalize()
def get_center_si_ti(self):
face, i, j, orientation = self.to_face_ij_orientation()
if self.is_leaf():
delta = 1
elif ((i ^ (self.id() >> 2)) & 1) != 0:
delta = 2
else:
delta = 0
return face, 2 * i + delta, 2 * j + delta
def get_center_uv(self):
face, si, ti = self.get_center_si_ti()
cls = self.__class__
return (cls.st_to_uv((0.5 / cls.MAX_SIZE) * si),
cls.st_to_uv((0.5 / cls.MAX_SIZE) * ti))
def to_face_ij_orientation(self):
i, j = 0, 0
face = self.face()
bits = (face & SWAP_MASK)
for k in xrange(7, -1, -1):
if k == 7:
nbits = (self.__class__.MAX_LEVEL - 7 * LOOKUP_BITS)
else:
nbits = LOOKUP_BITS
bits += ((self.id() >> (k * 2 * LOOKUP_BITS + 1) \
& ((1 << (2 * nbits)) - 1))) << 2
bits = LOOKUP_IJ[bits]
i += (bits >> (LOOKUP_BITS + 2)) << (k * LOOKUP_BITS)
j += ((bits >> 2) & ((1 << LOOKUP_BITS) - 1)) << (k * LOOKUP_BITS)
bits &= (SWAP_MASK | INVERT_MASK)
assert 0 == POS_TO_ORIENTATION[2]
assert SWAP_MASK == POS_TO_ORIENTATION[0]
if (self.lsb() & 0x1111111111111110) != 0:
bits ^= SWAP_MASK
orientation = bits
return face, i, j, orientation
def get_edge_neighbors(self):
level = self.level()
size = self.get_size_ij(level)
face, i, j, orientation = self.to_face_ij_orientation()
return (self.from_face_ij_same(
face, i, j - size, j - size >= 0).parent(level),
self.from_face_ij_same(
face, i + size, j,
i + size < self.__class__.MAX_SIZE).parent(level),
self.from_face_ij_same(
face, i, j + size,
j + size < self.__class__.MAX_SIZE).parent(level),
self.from_face_ij_same(
face, i - size, j, i - size >= 0).parent(level))
''' Return the neighbors of closest vertex to this cell at the given level.
Normally there are four neighbors, but the closest vertex may only have
three neighbors if it is one of the 8 cube vertices.
'''
def get_vertex_neighbors(self, level):
# "level" must be strictly less than this cell's level so that we can
# determine which vertex this cell is closest to.
assert level < self.level()
face, i, j, orientation = self.to_face_ij_orientation()
# Determine the i- and j-offsets to the closest neighboring cell in each
# direction. This involves looking at the next bit of "i" and "j" to
# determine which quadrant of this->parent(level) this cell lies in.
halfsize = self.get_size_ij(level + 1)
size = halfsize << 1
if i & halfsize:
ioffset = size
isame = (i + size) < self.__class__.MAX_SIZE
else:
ioffset = -size
isame = (i - size) >= 0
if j & halfsize:
joffset = size
jsame = (j + size) < self.__class__.MAX_SIZE
else:
joffset = -size
jsame = (j - size) >= 0
neighbors = []
neighbors.append(self.parent(level))
neighbors.append(self.__class__.from_face_ij_same(face, i + ioffset,
j, isame).parent(level))
neighbors.append(self.__class__.from_face_ij_same(face, i,
j + joffset, jsame).parent(level))
if isame or jsame:
neighbors.append(self.__class__.from_face_ij_same(face, i + ioffset,
j + joffset, isame and jsame).parent(level))
return neighbors
def get_all_neighbors(self, nbr_level):
face, i, j, orientation = self.to_face_ij_orientation()
# Find the coordinates of the lower left-hand leaf cell. We need to
# normalize (i,j) to a known position within the cell because nbr_level
# may be larger than this cell's level.
size = self.get_size_ij()
i &= -size
j &= -size
nbr_size = self.get_size_ij(nbr_level)
assert nbr_size <= size
# We compute the N-S, E-W, and diagonal neighbors in one pass.
# The loop test is at the end of the loop to avoid 32-bit overflow.
k = -nbr_size
while True:
if k < 0:
same_face = (j + k >= 0)
elif k >= size:
same_face = (j + k < self.__class__.MAX_SIZE)
else:
same_face = False
# North and South neighbors.
yield self.__class__.from_face_ij_same(
face, i + k, j - nbr_size,
j - size >= 0).parent(nbr_level)
yield self.__class__.from_face_ij_same(face, i + k, j + size,
j + size < self.__class__.MAX_SIZE).parent(nbr_level)
yield self.__class__.from_face_ij_same(face, i - nbr_size, j + k,
same_face and i - size >= 0).parent(nbr_level)
yield self.__class__.from_face_ij_same(face, i + size, j + k,
same_face and i + size < self.__class__.MAX_SIZE) \
.parent(nbr_level)
if k >= size:
break
k += nbr_size
def get_size_ij(self, *args):
if len(args) == 0:
level = self.level()
else:
level = args[0]
return 1 << (self.__class__.MAX_LEVEL - level)
def to_token(self):
return format(self.id(), 'x')
@classmethod
def from_token(cls, token):
return cls(int(token, 16))
@classmethod
def st_to_uv(cls, s):
if cls.PROJECTION == cls.LINEAR_PROJECTION:
return 2 * s - 1
elif cls.PROJECTION == cls.TAN_PROJECTION:
s = math.tan((math.pi / 2.0) * s - math.pi / 4.0)
return s + (1.0 / (1 << 53)) * s
elif cls.PROJECTION == cls.QUADRATIC_PROJECTION:
if s >= 0.5:
return (1.0 / 3.0) * (4 * s * s - 1)
else:
return (1.0 / 3.0) * (1 - 4 * (1 - s) * (1 - s))
else:
raise ValueError('unknown projection type')
@classmethod
def uv_to_st(cls, u):
if cls.PROJECTION == cls.LINEAR_PROJECTION:
return 0.5 * (u + 1)
elif cls.PROJECTION == cls.TAN_PROJECTION:
return (2 * (1.0 / math.pi)) * (math.atan(u) * math.pi / 4.0)
elif cls.PROJECTION == cls.QUADRATIC_PROJECTION:
if u >= 0:
return 0.5 * math.sqrt(1 + 3 * u)
else:
return 1 - 0.5 * math.sqrt(1 - 3 * u)
else:
raise ValueError('unknown projection type')
@classmethod
def avg_area(cls):
# Same for all projections.
return AreaMetric(4 * math.pi / 6)
@classmethod
def max_edge(cls):
return LengthMetric(cls.max_angle_span().deriv())
@classmethod
def max_angle_span(cls):
if cls.PROJECTION == cls.LINEAR_PROJECTION:
return LengthMetric(2)
elif cls.PROJECTION == cls.TAN_PROJECTION:
return LengthMetric(math.pi / 2)
elif cls.PROJECTION == cls.QUADRATIC_PROJECTION:
return LengthMetric(1.704897179199218452)
else:
raise ValueError('unknown projection type')
@classmethod
def max_diag(cls):
if cls.PROJECTION == cls.LINEAR_PROJECTION:
return LengthMetric(2 * math.sqrt(2))
elif cls.PROJECTION == cls.TAN_PROJECTION:
return LengthMetric(math.pi * math.sqrt(2.0 / 3.0))
elif cls.PROJECTION == cls.QUADRATIC_PROJECTION:
return LengthMetric(2.438654594434021032)
else:
raise ValueError('unknown projection type')
@classmethod
def min_width(cls):
if cls.PROJECTION == cls.LINEAR_PROJECTION:
return LengthMetric(math.sqrt(2))
elif cls.PROJECTION == cls.TAN_PROJECTION:
return LengthMetric(math.pi / 2 * math.sqrt(2))
elif cls.PROJECTION == cls.QUADRATIC_PROJECTION:
return LengthMetric(2 * math.sqrt(2) / 3)
else:
raise ValueError('unknown projection type')
class Metric(object):
def __init__(self, deriv, dim):
self.__deriv = deriv
self.__dim = dim
def deriv(self):
return self.__deriv
def get_value(self, level):
return math.ldexp(self.deriv(), -self.__dim * level)
def get_closest_level(self, value):
raise NotImplementedError()
def get_min_level(self, value):
raise NotImplementedError()
def get_max_level(self, value):
if value <= 0:
return CellId.MAX_LEVEL
m, x = math.frexp(self.deriv() / value)
level = max(0, min(CellId.MAX_LEVEL, (x - 1) >> (self.__dim - 1)))
assert level == 0 or self.get_value(level) >= value
assert level == CellId.MAX_LEVEL or self.get_value(level + 1) < value
return level
class LengthMetric(Metric):
def __init__(self, deriv):
super(LengthMetric, self).__init__(deriv, 1)
class AreaMetric(Metric):
def __init__(self, deriv):
super(AreaMetric, self).__init__(deriv, 2)
# like fmod but rounds to nearest integer instead of floor
def drem(x, y):
xd = decimal.Decimal(x)
yd = decimal.Decimal(y)
return float(xd.remainder_near(yd))
# functions originally in S2 source file
def valid_face_xyz_to_uv(face, p):
assert p.dot_prod(face_uv_to_xyz(face, 0, 0)) > 0
if face == 0:
return p[1] / p[0], p[2] / p[0]
elif face == 1:
return -p[0] / p[1], p[2] / p[1]
elif face == 2:
return -p[0] / p[2], -p[1] / p[2]
elif face == 3:
return p[2] / p[0], p[1] / p[0]
elif face == 4:
return p[2] / p[1], -p[0] / p[1]
else:
return -p[1] / p[2], -p[0] / p[2]
def xyz_to_face_uv(p):
face = p.largest_abs_component()
if p[face] < 0:
face += 3
u, v = valid_face_xyz_to_uv(face, p)
return face, u, v
def face_xyz_to_uv(face, p):
if face < 3:
if p[face] <= 0:
return False, 0, 0
else:
if p[face - 3] >= 0:
return False, 0, 0
u, v = valid_face_xyz_to_uv(face, p)
return True, u, v
def face_uv_to_xyz(face, u, v):
if face == 0:
return Point(1, u, v)
elif face == 1:
return Point(-u, 1, v)
elif face == 2:
return Point(-u, -v, 1)
elif face == 3:
return Point(-1, -v, -u)
elif face == 4:
return Point(v, -1, -u)
else:
return Point(v, u, -1)
def get_norm(face):
return face_uv_to_xyz(face, 0, 0)
# Return the right-handed normal (not necessarily unit length) for an
# edge in the direction of the positive v-axis at the given u-value on
# the given face. (This vector is perpendicular to the plane through
# the sphere origin that contains the given edge.)
def get_u_norm(face, u):
if face == 0:
return Point(u, -1, 0)
elif face == 1:
return Point(1, u, 0)
elif face == 2:
return Point(1, 0, u)
elif face == 3:
return Point(-u, 0, 1)
elif face == 4:
return Point(0, -u, 1)
else:
return Point(0, -1, -u)
# Return the right-handed normal (not necessarily unit length) for an
# edge in the direction of the positive u-axis at the given v-value on
# the given face.
def get_v_norm(face, v):
if face == 0:
return Point(-v, 0, 1)
elif face == 1:
return Point(0, -v, 1)
elif face == 2:
return Point(0, -1, -v)
elif face == 3:
return Point(v, -1, 0)
elif face == 4:
return Point(1, v, 0)
else:
return Point(1, 0, v)
def get_u_axis(face):
if face == 0:
return Point(0, 1, 0)
elif face == 1:
return Point(-1, 0, 0)
elif face == 2:
return Point(-1, 0, 0)
elif face == 3:
return Point(0, 0, -1)
elif face == 4:
return Point(0, 0, -1)
else:
return Point(0, 1, 0)
def get_v_axis(face):
if face == 0:
return Point(0, 0, 1)
elif face == 1:
return Point(0, 0, 1)
elif face == 2:
return Point(0, -1, 0)
elif face == 3:
return Point(0, -1, 0)
elif face == 4:
return Point(1, 0, 0)
else:
return Point(1, 0, 0)
def is_unit_length(p):
return math.fabs(p.norm() * p.norm() - 1) <= 1e-15
def ortho(a):
k = a.largest_abs_component() - 1
if k < 0:
k = 2
if k == 0:
temp = Point(1, 0.0053, 0.00457)
elif k == 1:
temp = Point(0.012, 1, 0.00457)
else:
temp = Point(0.012, 0.0053, 1)
return a.cross_prod(temp).normalize()
def origin():
# These values are ones that try not to overlap cells etc
return Point(0.00457, 1, 0.0321).normalize()
'''
The direction of a.CrossProd(b) becomes unstable as (a + b) or (a - b)
approaches zero. This leads to situations where a.CrossProd(b) is not
very orthogonal to "a" and/or "b". We could fix this using Gram-Schmidt,
but we also want b.RobustCrossProd(a) == -a.RobustCrossProd(b).
The easiest fix is to just compute the cross product of (b+a) and (b-a).
Mathematically, this cross product is exactly twice the cross product of
"a" and "b", but it has the numerical advantage that (b+a) and (b-a)
are always perpendicular (since "a" and "b" are unit length). This
yields a result that is nearly orthogonal to both "a" and "b" even if
these two values differ only in the lowest bit of one component.
'''
def robust_cross_prod(a, b):
assert is_unit_length(a)
assert is_unit_length(b)
x = (b + a).cross_prod(b - a)
if x != Point(0, 0, 0):
return x
return ortho(a)
def simple_crossing(a, b, c, d):
ab = a.cross_prod(b)
acb = -(ab.dot_prod(c))
bda = ab.dot_prod(d)
if acb * bda <= 0:
return False
cd = c.cross_prod(d)
cbd = -(cd.dot_prod(b))
dac = cd.dot_prod(a)
return (acb * cbd > 0) and (acb * dac > 0)
def girard_area(a, b, c):
ab = robust_cross_prod(a, b)
bc = robust_cross_prod(b, c)
ac = robust_cross_prod(a, c)
return max(0.0, ab.angle(ac) - ab.angle(bc) + bc.angle(ac))
def area(a, b, c):
assert is_unit_length(a)
assert is_unit_length(b)
assert is_unit_length(c)
sa = b.angle(c)
sb = c.angle(a)
sc = a.angle(b)
s = 0.5 * (sa + sb + sc)
if s >= 3e-4:
s2 = s * s
dmin = s - max(sa, max(sb, sc))
if dmin < 1e-2 * s * s2 * s2:
area = girard_area(a, b, c)
if dmin < 2 * (0.1 * area):
return area
return 4 * math.atan(math.sqrt(max(0.0, math.tan(0.5 * s) \
* math.tan(0.5 * (s - sa)) \
* math.tan(0.5 * (s - sb)) \
* math.tan(0.5 * (s - sc)))))
# Return true if the points A, B, C are strictly counterclockwise. Return
# false if the points are clockwise or collinear (i.e. if they are all
# contained on some great circle).
#
# Due to numerical errors, situations may arise that are mathematically
# impossible, e.g. ABC may be considered strictly CCW while BCA is not.
# However, the implementation guarantees the following:
#
# If simple_ccw(a,b,c), then !simple_ccw(c,b,a) for all a,b,c.
def simple_ccw(a, b, c):
# We compute the signed volume of the parallelepiped ABC. The usual
# formula for this is (AxB).C, but we compute it here using (CxA).B
# in order to ensure that ABC and CBA are not both CCW. This follows
# from the following identities (which are true numerically, not just
# mathematically):
#
# (1) x.CrossProd(y) == -(y.CrossProd(x))
# (2) (-x).DotProd(y) == -(x.DotProd(y))
return c.cross_prod(a).dot_prod(b) > 0
class Interval(object):
def __init__(self, lo, hi):
#self.__bounds = [args[0], args[1]]
self.__bounds = (lo, hi)
def __repr__(self):
return '{}: ({}, {})'.format(
self.__class__.__name__, self.__bounds[0], self.__bounds[1])
def lo(self):
return self.__bounds[0]
def hi(self):
return self.__bounds[1]
def bound(self, i):
return self.__bounds[i]
def bounds(self):
return self.__bounds
@classmethod
def empty(cls):
return cls()
class LineInterval(Interval):
def __init__(self, lo=1, hi=0):
super(LineInterval, self).__init__(lo, hi)
def __eq__(self, other):
return isinstance(other, self.__class__) \
and ((self.lo() == other.lo() and self.hi() == other.hi()) \
or (self.is_empty() and other.is_empty()))
def __hash__(self):
return hash(self.__bounds)
@classmethod
def from_point_pair(cls, a, b):
if a <= b:
return cls(a, b)
else:
return cls(b, a)
def contains(self, other):
if isinstance(other, self.__class__):
if other.is_empty():
return True
return other.lo() >= self.lo() and other.hi() <= self.hi()
else:
return other >= self.lo() and other <= self.hi()
def interior_contains(self, other):
if isinstance(other, self.__class__):
if other.is_empty():
return True
return other.lo() > self.lo() and other.hi() < self.hi()
else:
return other > self.lo() and other < self.hi()
def intersects(self, other):
if self.lo() <= other.lo():
return other.lo() <= self.hi() and other.lo() <= other.hi()
else:
return self.lo() <= other.hi() and self.lo() <= self.hi()
def interior_intersects(self, other):
return other.lo() < self.hi() and self.lo() < other.hi() \
and self.lo() < self.hi() and other.lo() <= other.hi()
def union(self, other):
if self.is_empty():
return other
if other.is_empty():
return self
return self.__class__(min(self.lo(), other.lo()),
max(self.hi(), other.hi()))
def intersection(self, other):
return self.__class__(max(self.lo(), other.lo()),
min(self.hi(), other.hi()))
def expanded(self, radius):
assert radius >= 0
if self.is_empty():
return self
return self.__class__(self.lo() - radius, self.hi() + radius)
def get_center(self):
return 0.5 * (self.lo() + self.hi())
def get_length(self):
return self.hi() - self.lo()
def is_empty(self):
return self.lo() > self.hi()
def approx_equals(self, other, max_error=1e-15):
if self.is_empty():
return other.get_length() <= max_error
if other.is_empty():
return self.get_length() <= max_error
return math.fabs(other.lo() - self.lo()) \
+ math.fabs(other.hi() - self.hi()) <= max_error
# originally S1Interval in C++ code
class SphereInterval(Interval):
def __init__(self, lo=math.pi, hi=-math.pi, args_checked=False):
if args_checked:
super(SphereInterval, self).__init__(lo, hi)
else:
clamped_lo, clamped_hi = lo, hi
if lo == -math.pi and hi != math.pi:
clamped_lo = math.pi
if hi == -math.pi and lo != math.pi:
clamped_hi = math.pi
super(SphereInterval, self).__init__(clamped_lo, clamped_hi)
assert self.is_valid()
def __eq__(self, other):
return isinstance(other, self.__class__) \
and self.lo() == other.lo() and self.hi() == other.hi()
@classmethod
def from_point_pair(cls, a, b):
assert math.fabs(a) <= math.pi
assert math.fabs(b) <= math.pi
if a == -math.pi:
a = math.pi
if b == -math.pi:
b = math.pi
if cls.positive_distance(a, b) <= math.pi:
return cls(a, b, args_checked=True)
else:
return cls(b, a, args_checked=True)
@classmethod
def positive_distance(cls, a, b):
d = b - a
if d >= 0:
return d
return (b + math.pi) - (a - math.pi)
@classmethod
def full(cls):
return cls(-math.pi, math.pi, args_checked=True)
def is_full(self):
return (self.hi() - self.lo()) == (2 * math.pi)
def is_valid(self):
return math.fabs(self.lo()) <= math.pi \
and math.fabs(self.hi()) <= math.pi \
and not (self.lo() == -math.pi and self.hi() != math.pi) \
and not (self.hi() == -math.pi and self.lo() != math.pi)
def is_inverted(self):
return self.lo() > self.hi()
def is_empty(self):
return self.lo() - self.hi() == 2 * math.pi
def get_center(self):
center = 0.5 * (self.lo() + self.hi())
if not self.is_inverted():
return center
if center <= 0:
return center + math.pi
else:
return center - math.pi
def get_length(self):
length = self.hi() - self.lo()
if length >= 0:
return length
length += (2 * math.pi)
if length > 0:
return length
else:
return -1
# Return the complement of the interior of the interval. An interval and
# its complement have the same boundary but do not share any interior
# values. The complement operator is not a bijection, since the complement
# of a singleton interval (containing a single value) is the same as the
# complement of an empty interval.
def complement(self):
if self.lo() == self.hi():
return self.__class__.full()
return self.__class__(self.hi(), self.lo())
def approx_equals(self, other, max_error=1e-15):
if self.is_empty():
return other.get_length() <= max_error
if other.is_empty():
return self.get_length() <= max_error
return (math.fabs(drem(other.lo() - self.lo(), 2 * math.pi)) \
+ math.fabs(drem(other.hi() - self.hi(), 2 * math.pi))) <= max_error
def fast_contains(self, other):
if self.is_inverted():
return (other >= self.lo() or other <= self.hi()) \
and not self.is_empty()
else:
return other >= self.lo() and other <= self.hi()
def contains(self, other):
if isinstance(other, self.__class__):
if self.is_inverted():
if other.is_inverted():
return other.lo() >= self.lo() and other.hi() <= self.hi()
return (other.lo() >= self.lo() or other.hi() <= self.hi()) \
and not self.is_empty()
else:
if other.is_inverted():
return self.is_full() or other.is_empty()
return other.lo() >= self.lo() and other.hi() <= self.hi()
else:
assert math.fabs(other) <= math.pi
if other == -math.pi:
other = math.pi
return self.fast_contains(other)
def interior_contains(self, other):
if isinstance(other, self.__class__):
if self.is_inverted():
if not other.is_inverted():
return other.lo() > self.lo() or other.hi() < self.hi()
return (other.lo() > self.lo() and other.hi() < self.hi()) \
or other.is_empty()
else:
if other.is_inverted():
return self.is_full() or other.is_empty()
return (other.lo() > self.lo() and other.hi() < self.hi()) \
or self.is_full()
else:
assert math.fabs(other) <= math.pi
if other == -math.pi:
other = math.pi
if self.is_inverted():
return other > self.lo() or other < self.hi()
else:
return (other > self.lo() and other < self.hi()) \
or self.is_full()
def intersects(self, other):
if self.is_empty() or other.is_empty():
return False
if self.is_inverted():
return other.is_inverted() or other.lo() <= self.hi() \
or other.hi() >= self.lo()
else:
if other.is_inverted():
return other.lo() <= self.hi() or other.hi() >= self.lo()
return other.lo() <= self.hi() and other.hi() >= self.lo()
def interior_intersects(self, other):
if self.is_empty() or other.is_empty() or self.lo() == self.hi():
return False
if self.is_inverted():
return other.is_inverted() or other.lo() < self.hi() \
or other.hi() > self.lo()
else:
if other.is_inverted():
return other.lo() < self.hi() or other.hi() > self.lo()
return (other.lo() < self.hi() and other.hi() > self.lo()) \
or self.is_full()
def union(self, other):
if other.is_empty():
return self
if self.fast_contains(other.lo()):
if self.fast_contains(other.hi()):
if self.contains(other):
return self
return self.__class__.full()
return self.__class__(self.lo(), other.hi(), args_checked=True)
if self.fast_contains(other.hi()):
return self.__class__(other.lo(), self.hi(), args_checked=True)
if self.is_empty() or other.fast_contains(self.lo()):
return other
dlo = self.__class__.positive_distance(other.hi(), self.lo())
dhi = self.__class__.positive_distance(self.hi(), other.lo())
if dlo < dhi:
return self.__class__(other.lo(), self.hi(), args_checked=True)
else:
return self.__class__(self.lo(), other.hi(), args_checked=True)
def intersection(self, other):
if other.is_empty():
return self.__class__.empty()
if self.fast_contains(other.lo()):
if self.fast_contains(other.hi()):
if other.get_length() < self.get_length():
return other
return self
return self.__class__(other.lo(), self.hi(), args_checked=True)
if self.fast_contains(other.hi()):
return self.__class__(self.lo(), other.hi(), args_checked=True)
if other.fast_contains(self.lo()):
return self
assert not self.intersects(other)
return self.__class__.empty()
def expanded(self, radius):
assert radius >= 0
if self.is_empty():
return self
if self.get_length() + 2 * radius >= 2 * math.pi - 1e-15:
return self.__class__.full()
lo = drem(self.lo() - radius, 2 * math.pi)
hi = drem(self.hi() + radius, 2 * math.pi)
if lo <= -math.pi:
lo = math.pi
return self.__class__(lo, hi)
def get_complement_center(self):
if self.lo() != self.hi():
return self.complement().get_center()
else:
if self.hi() <= 0:
return self.hi() + math.pi
else:
return self.hi() - math.pi
def get_directed_hausdorff_distance(self, other):
if other.contains(self):
return 0.0
if other.is_empty():
return math.pi
other_complement_center = other.get_complement_center()
if self.contains(other_complement_center):
return self.__class__.positive_distance(other.hi(),
other_complement_center)
else:
if self.__class__(other.hi(), other_complement_center) \
.contains(self.hi()):
hi_hi = self.__class__.positive_distance(other.hi(), self.hi())
else:
hi_hi = 0
if self.__class__(other_complement_center, other.lo()) \
.contains(self.lo()):
lo_lo = self.__class__.positive_distance(self.lo(), other.lo())
else:
lo_lo = 0
assert hi_hi > 0 or lo_lo > 0
return max(hi_hi, lo_lo)
class Cell(object):
def __init__(self, cell_id=None):
self.__uv = [[None, None], [None, None]]
if cell_id is not None:
self.__cell_id = cell_id
face, i, j, orientation = cell_id.to_face_ij_orientation()
ij = (i, j)
self.__face = face
self.__orientation = orientation
self.__level = cell_id.level()
cell_size = cell_id.get_size_ij()
for d in xrange(2):
ij_lo = ij[d] & -cell_size
ij_hi = ij_lo + cell_size
self.__uv[d][0] = CellId.st_to_uv(
(1.0 / CellId.MAX_SIZE) * ij_lo)
self.__uv[d][1] = CellId.st_to_uv(
(1.0 / CellId.MAX_SIZE) * ij_hi)
@classmethod
def from_lat_lon(cls, lat_lon):
return cls(CellId.from_lat_lon(lat_lon))
@classmethod
def from_point(cls, point):
return cls(CellId.from_point(point))
def __repr__(self):
return '{}: face {}, level {}, orientation {}, id {}'.format(
self.__class__.__name__, self.__face, self.__level,
self.orientation(), self.__cell_id.id())
@classmethod
def from_face_pos_level(cls, face, pos, level):
return cls(CellId.from_face_pos_level(face, pos, level))
def id(self):
return self.__cell_id
def face(self):
return self.__face
def level(self):
return self.__level
def orientation(self):
return self.__orientation
def is_leaf(self):
return self.__level == CellId.MAX_LEVEL
# Return the inward-facing normal of the great circle passing through
# the edge from vertex k to vertex k+1 (mod 4). The normals returned
# by GetEdgeRaw are not necessarily unit length.
def get_edge(self, k):
return self.get_edge_raw(k).normalize()
def get_edge_raw(self, k):
if k == 0:
return get_v_norm(self.__face, self.__uv[1][0]) # South
elif k == 1:
return get_u_norm(self.__face, self.__uv[0][1]) # East
elif k == 2:
return -get_v_norm(self.__face, self.__uv[1][1]) # North
else:
return -get_u_norm(self.__face, self.__uv[0][0]) # West
# Return the k-th vertex of the cell (k = 0,1,2,3). Vertices are returned
# in CCW order. The points returned by GetVertexRaw are not necessarily
# unit length.
def get_vertex(self, k):
return self.get_vertex_raw(k).normalize()
def get_vertex_raw(self, k):
return face_uv_to_xyz(self.__face,
self.__uv[0][(k >> 1) ^ (k & 1)], self.__uv[1][k >> 1])
# Return the area of this cell as accurately as possible. This method is
# more expensive but it is accurate to 6 digits of precision even for leaf
# cells (whose area is approximately 1e-18).
def exact_area(self):
v0 = self.get_vertex(0)
v1 = self.get_vertex(1)
v2 = self.get_vertex(2)
v3 = self.get_vertex(3)
return area(v0, v1, v2) + area(v0, v2, v3)
def average_area(self):
return CellId.avg_area().get_value(self.__level)
def approx_area(self):
if self.__level < 2:
return self.average_area()
flat_area = 0.5 * (self.get_vertex(2) - self.get_vertex(0)) \
.cross_prod(self.get_vertex(3) - self.get_vertex(1)) \
.norm()
return flat_area * 2 / (1 + math.sqrt(1 - min(
(1.0 / math.pi) * flat_area, 1.0)))
def subdivide(self):
uv_mid = self.__cell_id.get_center_uv()
for pos, cell_id in enumerate(self.__cell_id.children()):
child = Cell()
child.__face = self.__face
child.__level = self.__level + 1
child.__orientation = self.__orientation ^ POS_TO_ORIENTATION[pos]
child.__cell_id = cell_id
# We want to split the cell in half in "u" and "v". To decide which
# side to set equal to the midpoint value, we look at cell's (i,j)
# position within its parent. The index for "i" is in bit 1 of ij.
ij = POS_TO_IJ[self.__orientation][pos]
i = ij >> 1
j = ij & 1
child.__uv[0][i] = self.__uv[0][i]
child.__uv[0][1 - i] = uv_mid[0]
child.__uv[1][j] = self.__uv[1][j]
child.__uv[1][1 - j] = uv_mid[1]
yield child
def get_center(self):
return self.get_center_raw().normalize()
def get_center_raw(self):
return self.__cell_id.to_point_raw()
def contains(self, other):
if isinstance(other, self.__class__):
return self.__cell_id.contains(other.__cell_id)
elif isinstance(other, Point):
valid, u, v = face_xyz_to_uv(self.__face, other)
if not valid:
return False
return u >= self.__uv[0][0] and u <= self.__uv[0][1] \
and v >= self.__uv[1][0] and v <= self.__uv[1][1]
def may_intersect(self, cell):
return self.__cell_id.intersects(cell.__cell_id)
def get_latitude(self, i, j):
p = face_uv_to_xyz(self.__face, self.__uv[0][i], self.__uv[1][j])
return LatLon.latitude(p).radians
def get_longitude(self, i, j):
p = face_uv_to_xyz(self.__face, self.__uv[0][i], self.__uv[1][j])
return LatLon.longitude(p).radians
def get_cap_bound(self):
u = 0.5 * (self.__uv[0][0] + self.__uv[0][1])
v = 0.5 * (self.__uv[1][0] + self.__uv[1][1])
cap = Cap.from_axis_height(
face_uv_to_xyz(self.__face, u, v).normalize(), 0)
for k in xrange(4):
cap.add_point(self.get_vertex(k))
return cap
def get_rect_bound(self):
if self.__level > 0:
u = self.__uv[0][0] + self.__uv[0][1]
v = self.__uv[1][0] + self.__uv[1][1]
if get_u_axis(self.__face)[2] == 0:
i = int(u < 0)
else:
i = int(u > 0)
if get_v_axis(self.__face)[2] == 0:
j = int(v < 0)
else:
j = int(v > 0)
max_error = 1.0 / (1 << 51)
lat = LineInterval.from_point_pair(self.get_latitude(i, j),
self.get_latitude(1 - i, 1 - j))
lat = lat.expanded(max_error).intersection(LatLonRect.full_lat())
if lat.lo() == (-math.pi / 2.0) or lat.hi() == (math.pi / 2.0):
return LatLonRect(lat, SphereInterval.full())
lon = SphereInterval.from_point_pair(self.get_longitude(i, 1 - j),
self.get_longitude(1 - i, j))
return LatLonRect(lat, lon.expanded(max_error))
pole_min_lat = math.asin(math.sqrt(1.0 / 3.0))
if self.__face == 0:
return LatLonRect(
LineInterval(-math.pi / 4.0, math.pi / 4.0),
SphereInterval(-math.pi / 4.0, math.pi / 4.0))
elif self.__face == 1:
return LatLonRect(
LineInterval(-math.pi / 4.0, math.pi / 4.0),
SphereInterval(math.pi / 4.0, 3.0 * math.pi / 4.0))
elif self.__face == 2:
return LatLonRect(
LineInterval(pole_min_lat, math.pi / 2.0),
SphereInterval(-math.pi, math.pi))
elif self.__face == 3:
return LatLonRect(
LineInterval(-math.pi / 4.0, math.pi / 4.0),
SphereInterval(3.0 * math.pi / 4.0, -3.0 * math.pi / 4.0))
elif self.__face == 4:
return LatLonRect(
LineInterval(-math.pi / 4.0, math.pi / 4.0),
SphereInterval(-3.0 * math.pi / 4.0, -math.pi / 4.0))
else:
return LatLonRect(
LineInterval(-math.pi / 2.0, -pole_min_lat),
SphereInterval(-math.pi, math.pi))
class CellUnion(object):
def __init__(self, cell_ids=None, raw=True):
if cell_ids is None:
self.__cell_ids = []
else:
self.__cell_ids = [cell_id
if isinstance(cell_id, CellId)
else CellId(cell_id)
for cell_id in cell_ids]
if raw:
self.normalize()
def __eq__(self, other):
return isinstance(other, self.__class__) \
and self.__cell_ids == other.__cell_ids
def __hash__(self):
return hash(self.__cell_ids)
def __repr__(self):
return '{}: {}'.format(self.__class__.__name__, self.__cell_ids)
@classmethod
def get_union(cls, x, y):
return cls(x.__cell_ids + y.__cell_ids)
@classmethod
def get_intersection(cls, *args):
if len(args) == 2 and isinstance(args[0], cls) \
and isinstance(args[1], CellId):
cell_union, cell_id = args
if cell_union.contains(cell_id):
return cls([cell_id])
else:
index = bisect.bisect_left(cell_union.__cell_ids,
cell_id.range_min())
idmax = cell_id.range_max()
intersected_cell_ids = []
while index != len(cell_union.__cell_ids) \
and cell_union.__cell_ids[index] <= idmax:
intersected_cell_ids.append(cell_union.__cell_ids[index])
index += 1
return cls(intersected_cell_ids)
elif len(args) == 2 and isinstance(args[0], cls) \
and isinstance(args[1], cls):
x, y = args
i, j = 0, 0
cell_ids = []
while i < x.num_cells() and j < y.num_cells():
imin = x.__cell_ids[i].range_min()
jmin = y.__cell_ids[j].range_min()
if imin > jmin:
if x.__cell_ids[i] <= y.__cell_ids[j].range_max():
cell_ids.append(x.__cell_ids[i])
i += 1
else:
j = bisect.bisect_left(y.__cell_ids, imin, lo=j + 1)
if x.__cell_ids[i] <= y.__cell_ids[j - 1].range_max():
j -= 1
elif jmin > imin:
if y.__cell_ids[j] <= x.__cell_ids[i].range_max():
cell_ids.append(y.__cell_ids[j])
j += 1
else:
i = bisect.bisect_left(x.__cell_ids, jmin, lo=i + 1)
if y.__cell_ids[j] <= x.__cell_ids[i - 1].range_max():
i -= 1
else:
if x.__cell_ids[i] < y.__cell_ids[j]:
cell_ids.append(x.__cell_ids[i])
i += 1
else:
cell_ids.append(y.__cell_ids[j])
j += 1
cell_union = cls(cell_ids)
assert all(cell_union.__cell_ids[i] <= cell_union.__cell_ids[i + 1]
for i in xrange(cell_union.num_cells() - 1))
assert not cell_union.normalize()
return cell_union
else:
raise NotImplementedError()
def expand(self, *args):
if len(args) == 1 and isinstance(args[0], (int, long)):
level = args[0]
output = []
level_lsb = CellId.lsb_for_level(level)
i = self.num_cells() - 1
while i >= 0:
cell_id = self.__cell_ids[i]
if cell_id.lsb() < level_lsb:
cell_id = cell_id.parent(level)
while i > 0 and cell_id.contains(self.__cell_ids[i - 1]):
i -= 1
output.append(cell_id)
cell_id.append_all_neighbors(level, output)
i -= 1
self.__cell_ids = output
elif len(args) == 2 and isinstance(args[0], Angle) \
and isinstance(args[1], (int, long)):
min_radius, max_level_diff = args
min_level = CellId.MAX_LEVEL
for cell_id in self.__cell_ids:
min_level = min(min_level, cell_id.level())
radius_level = CellId.min_width().get_max_level(
min_radius.radians)
if radius_level == 0 \
and min_radius.radians > CellId.min_width().get_value(0):
self.expand(0)
self.expand(min(min_level + max_level_diff, radius_level))
else:
raise NotImplementedError()
@classmethod
def get_difference(cls, x, y):
cell_ids = []
for cell_id in x.__cell_ids:
cls.__get_difference(cell_id, y, cell_ids)
cell_union = cls(cell_ids)
assert all(cell_union.__cell_ids[i] <= cell_union.__cell_ids[i + 1]
for i in xrange(cell_union.num_cells() - 1))
assert not cell_union.normalize()
return cell_union
@classmethod
def __get_difference(cls, cell_id, y, cell_ids):
if not y.intersects(cell_id):
cell_ids.append(cell_id)
elif not y.contains(cell_id):
for child in cell_id.children():
cls.__get_difference(child, y, cell_ids)
def num_cells(self):
return len(self.__cell_ids)
def cell_id(self, i):
return self.__cell_ids[i]
def cell_ids(self):
return self.__cell_ids
def normalize(self):
self.__cell_ids.sort()
output = []
for cell_id in self.__cell_ids:
if output and output[-1].contains(cell_id):
continue
while output and cell_id.contains(output[-1]):
output.pop()
while len(output) >= 3:
if (output[-3].id() ^ output[-2].id() ^ output[-1].id()) \
!= cell_id.id():
break
mask = cell_id.lsb() << 1
mask = ~(mask + (mask << 1))
id_masked = (cell_id.id() & mask)
if (output[-3].id() & mask) != id_masked \
or (output[-2].id() & mask) != id_masked \
or (output[-1].id() & mask) != id_masked \
or cell_id.is_face():
break
output.pop()
output.pop()
output.pop()
cell_id = cell_id.parent()
output.append(cell_id)
if len(output) < self.num_cells():
self.__cell_ids = output
return True
return False
def denormalize(self, min_level, level_mod):
assert min_level >= 0
assert min_level <= CellId.MAX_LEVEL
assert level_mod >= 1
assert level_mod <= 3
cell_ids = []
for cell_id in self.__cell_ids:
level = cell_id.level()
new_level = max(min_level, level)
if level_mod > 1:
new_level += (CellId.MAX_LEVEL - (new_level - min_level)) \
% level_mod
new_level = min(CellId.MAX_LEVEL, new_level)
if new_level == level:
cell_ids.append(cell_id)
else:
end = cell_id.child_end(new_level)
for child in cell_id.children(new_level):
cell_ids.append(child)
return cell_ids
def contains(self, *args):
if len(args) == 1 and isinstance(args[0], Cell):
return self.contains(args[0].id())
elif len(args) == 1 and isinstance(args[0], CellId):
cell_id = args[0]
index = bisect.bisect_left(self.__cell_ids, cell_id)
if index < len(self.__cell_ids) \
and self.__cell_ids[index].range_min() <= cell_id:
return True
return index != 0 \
and self.__cell_ids[index - 1].range_max() >= cell_id
elif len(args) == 1 and isinstance(args[0], Point):
return self.contains(CellId.from_point(args[0]))
elif len(args) == 1 and isinstance(args[0], self.__class__):
cell_union = args[0]
for i in xrange(cell_union.num_cells()):
if not self.contains(cell_union.cell_id(i)):
return False
return True
else:
raise NotImplementedError()
def intersects(self, *args):
if len(args) == 1 and isinstance(args[0], CellId):
cell_id = args[0]
index = bisect.bisect_left(self.__cell_ids, cell_id)
if index != len(self.__cell_ids) \
and self.__cell_ids[index].range_min() <= cell_id.range_max():
return True
return index != 0 \
and self.__cell_ids[index - 1].range_max() \
>= cell_id.range_min()
elif len(args) == 1 and isinstance(args[0], CellUnion):
cell_union = args[0]
for cell_id in cell_union.__cell_ids:
if self.intersects(cell_id):
return True
return False
else:
raise NotImplementedError()
FACE_CELLS = (Cell.from_face_pos_level(0, 0, 0),
Cell.from_face_pos_level(1, 0, 0),
Cell.from_face_pos_level(2, 0, 0),
Cell.from_face_pos_level(3, 0, 0),
Cell.from_face_pos_level(4, 0, 0),
Cell.from_face_pos_level(5, 0, 0))
class RegionCoverer(object):
class Candidate(object):
@property
def num_children(self):
return len(self.children)
def __init__(self):
self.__min_level = 0
self.__max_level = CellId.MAX_LEVEL
self.__level_mod = 1
self.__max_cells = 8
self.__region = None
self.__result = []
self.__pq = []
@property
def min_level(self):
return self.__min_level
@min_level.setter
def min_level(self, value):
assert value >= 0
assert value <= CellId.MAX_LEVEL
self.__min_level = max(0, min(CellId.MAX_LEVEL, value))
@property
def max_level(self):
return self.__max_level
@max_level.setter
def max_level(self, value):
assert value >= 0
assert value <= CellId.MAX_LEVEL
self.__max_level = max(0, min(CellId.MAX_LEVEL, value))
@property
def level_mod(self):
return self.__level_mod
@level_mod.setter
def level_mod(self, value):
assert value >= 1
assert value <= 3
self.__level_mod = max(1, min(3, value))
@property
def max_cells(self):
return self.__max_cells
@max_cells.setter
def max_cells(self, value):
self.__max_cells = value
def get_covering(self, region):
self.__result = []
tmp_union = self.__get_cell_union(region)
return tmp_union.denormalize(self.__min_level, self.__level_mod)
def get_interior_covering(self, region):
self.__result = []
tmp_union = self.__get_interior_cell_union(region)
return tmp_union.denormalize(self.__min_level, self.__level_mod)
def __new_candidate(self, cell):
if not self.__region.may_intersect(cell):
return None
is_terminal = False
if cell.level() >= self.__min_level:
if self.__interior_covering:
if self.__region.contains(cell):
is_terminal = True
elif cell.level() + self.__level_mod > self.__max_level:
return None
else:
if cell.level() + self.__level_mod > self.__max_level \
or self.__region.contains(cell):
is_terminal = True
candidate = self.__class__.Candidate()
candidate.cell = cell
candidate.is_terminal = is_terminal
candidate.children = []
return candidate
def __max_children_shift(self):
return 2 * self.__level_mod
def __expand_children(self, candidate, cell, num_levels):
num_levels -= 1
num_terminals = 0
for child_cell in cell.subdivide():
if num_levels > 0:
if self.__region.may_intersect(child_cell):
num_terminals += self.__expand_children(candidate,
child_cell, num_levels)
continue
child = self.__new_candidate(child_cell)
if child is not None:
candidate.children.append(child)
if child.is_terminal:
num_terminals += 1
return num_terminals
def __add_candidate(self, candidate):
if candidate is None:
return
if candidate.is_terminal:
self.__result.append(candidate.cell.id())
return
assert candidate.num_children == 0
num_levels = self.__level_mod
if candidate.cell.level() < self.__min_level:
num_levels = 1
num_terminals = self.__expand_children(candidate,
candidate.cell, num_levels)
if candidate.num_children == 0:
''' Not needed due to GC '''
elif not self.__interior_covering \
and num_terminals == 1 << self.__max_children_shift() \
and candidate.cell.level() >= self.__min_level:
candidate.is_terminal = True
self.__add_candidate(candidate)
else:
priority = (((
(candidate.cell.level() << self.__max_children_shift())
+ candidate.num_children) << self.__max_children_shift())
+ num_terminals)
heapq.heappush(self.__pq, (priority, candidate))
def __get_initial_candidates(self):
if self.__max_cells >= 4:
cap = self.__region.get_cap_bound()
level = min(CellId.min_width().get_max_level(
2 * cap.angle().radians),
min(self.__max_level, CellId.MAX_LEVEL - 1))
if self.__level_mod > 1 and level > self.__min_level:
level -= (level - self.__min_level) % self.__level_mod
if level > 0:
cell_id = CellId.from_point(cap.axis())
base = cell_id.get_vertex_neighbors(level)
for i in xrange(len(base)):
self.__add_candidate(self.__new_candidate(Cell(base[i])))
return
for face in xrange(6):
self.__add_candidate(self.__new_candidate(FACE_CELLS[face]))
def __get_covering(self, region):
assert len(self.__pq) == 0
assert len(self.__result) == 0
self.__region = region
self.__get_initial_candidates()
while len(self.__pq) > 0 \
and (not self.__interior_covering \
or len(self.__result) < self.__max_cells):
candidate = heapq.heappop(self.__pq)[1]
if self.__interior_covering:
result_size = 0
else:
result_size = len(self.__pq)
if candidate.cell.level() < self.__min_level \
or candidate.num_children == 1 \
or len(self.__result) + result_size + candidate.num_children \
<= self.__max_cells:
for child in candidate.children:
self.__add_candidate(child)
elif self.__interior_covering:
''' Do nothing here '''
else:
candidate.is_terminal = True
self.__add_candidate(candidate)
self.__pq = []
self.__region = None
def __get_cell_union(self, region):
self.__interior_covering = False
self.__get_covering(region)
return CellUnion(self.__result)
def __get_interior_cell_union(self, region):
self.__interior_covering = True
self.__get_covering(region)
return CellUnion(self.__result)
@classmethod
def flood_fill(cls, region, start):
all = set()
frontier = []
all.add(start)
frontier.append(start)
while len(frontier) != 0:
cell_id = frontier.pop()
if not region.may_intersect(Cell(cell_id)):
continue
yield cell_id
neighbors = cell_id.get_edge_neighbors()
for nbr in neighbors:
if nbr not in all:
all.add(nbr)
frontier.append(nbr)
@classmethod
def get_simple_covering(cls, region, start, level):
return cls.flood_fill(region, CellId.from_point(start).parent(level))
