EDIT3:
As it turns out, considering all possible variations is more complicated than I thought at first glance. This third iteration of my code should be correct for all possible inputs. Due to the increased complexity I have discarded the vectorized numpy variant. The generator version is the following:
def overlapping_sectors3(sectors, interval):
"""
Yields overlapping radial intervals.
Returns the overlapping intervals between each of the sector-intervals
and the comparison-interval.
Args:
sectors: List of intervals.
Interval borders must be in [0, 2*pi).
interval: Single interval aginst which the overlap is calculated.
Interval borders must be in [0, 2*pi).
Yields:
A list of intervals marking the overlaping areas.
Interval borders are guaranteed to be in [0, 2*pi).
"""
i_lhs, i_rhs = interval
if i_lhs > i_rhs:
for s_lhs, s_rhs in sectors:
if s_lhs > s_rhs:
# CASE 1
o_lhs = max(s_lhs, i_lhs)
# o_rhs = min(s_rhs+2*np.pi, i_rhs+2*np.pi)
o_rhs = min(s_rhs, i_rhs)
# since o_rhs > 2pi > o_lhs
yield o_lhs, o_rhs
#o_lhs = max(s_lhs+2pi, i_lhs)
# o_rhs = min(s_rhs+4pi, i_rhs+2pi)
# since o_lhs and o_rhs > 2pi
o_lhs = s_lhs
o_rhs = i_rhs
if o_lhs < o_rhs:
yield o_lhs, o_rhs
else:
# CASE 2
o_lhs = max(s_lhs, i_lhs)
# o_rhs = min(s_rhs, i_rhs+2*np.pi)
o_rhs = s_rhs # since i_rhs + 2pi > 2pi > s_rhs
if o_lhs < o_rhs:
yield o_lhs, o_rhs
# o_lhs = max(s_lhs+2pi, i_lhs)
# o_rhs = min(s_rhs+2pi, i_rhs+2pi)
# since s_lhs+2pi > 2pi > i_lhs and both o_lhs and o_rhs > 2pi
o_lhs = s_lhs
o_rhs = min(s_rhs, i_rhs)
if o_lhs < o_rhs:
yield o_lhs, o_rhs
else:
for s_lhs, s_rhs in sectors:
if s_lhs > s_rhs:
# CASE 3
o_lhs = max(s_lhs, i_lhs)
o_rhs = i_rhs
if o_lhs < o_rhs:
yield o_lhs, o_rhs
o_lhs = i_lhs
o_rhs = min(s_rhs, i_rhs)
if o_lhs < o_rhs:
yield o_lhs, o_rhs
else:
# CASE 4
o_lhs = max(s_lhs, i_lhs)
o_rhs = min(s_rhs, i_rhs)
if o_lhs < o_rhs:
yield o_lhs, o_rhs
It can be tested with:
import numpy as np
from collections import namedtuple
TestCase = namedtuple('TestCase', ['sectors', 'interval', 'expected', 'remark'])
testcases = []
def newcase(sectors, interval, expected, remark=None):
testcases.append( TestCase(sectors, interval, expected, remark) )
newcase(
[[280,70]],
[270,90],
[[280,70]],
"type 1"
)
newcase(
[[10,150]],
[270,90],
[[10,90]],
"type 2"
)
newcase(
[[10,150]],
[270,350],
[],
"type 4"
)
newcase(
[[50,350]],
[10,90],
[[50,90]],
"type 4"
)
newcase(
[[30,0]],
[300,60],
[[30,60],[300,0]],
"type 1"
)
newcase(
[[30,5]],
[300,60],
[[30,60],[300,5]],
"type 1"
)
newcase(
[[30,355]],
[300,60],
[[30,60],[300,355]],
"type 3"
)
def isequal(A,B):
if len(A) != len(B):
return False
A = np.array(A).round()
B = np.array(B).round()
a = set(map(tuple, A))
b = set(map(tuple, B))
return a == b
for caseindex, case in enumerate(testcases):
print("### testcase %2d ###" % caseindex)
print("sectors : %s" % case.sectors)
print("interval: %s" % case.interval)
if case.remark:
print(case.remark)
sectors = np.array(case.sectors)/180*np.pi
interval = np.array(case.interval)/180*np.pi
result = overlapping_sectors3(sectors, interval)
result = np.array(list(result))*180/np.pi
if isequal(case.expected, result):
print('PASS')
else:
print('FAIL')
print('\texp: %s' % case.expected)
print('\tgot: %s' % result)
To understand the logic behind it consider the following:
- Each interval has a left hand side (lhs) and a right hand side (rhs)
- if lhs > rhs then the interval "wraps round", i.e. it is actually the interval [lhs, rhs+2pi]
- when comparing the current sector aginst the comparison-interval we have to consider four cases
- both wrap round
- only the comparison-interval wraps round
- only the sector-interval wraps round
- neither one wraps round
- With ordinary intervals the overlapping interval is
[o_lhs, o_rhs] with o_lhs=max(lhs_1, lhs2) and o_rhs=min(rhs_1, rhs_2) iff o_lhs < o_rhs
- "Unwarpping" all intervals by adding
2pi to the rhs iff rhs<lhs yields intervals in [0, 4*np.pi)
- We will call
[0,2*pi) the first, [2*pi, 4*pi) the second and [4*pi, 6*pi) the third circumvolution.
The four cases:
- case 4: neither interval wraps round, so all borders lie within the first circumvolution. We can just calculate the overlap as with any orinariy intervals.
- case 2 and 3: exactly one interval wraps round. That means one interval (lets call it a) is entirely within the first circumvolution, while second one (lets call it b) spawns both the first and the second circumvolution. That means, that a can intersect b in both the first and the second circumvolution. First we consider the first circumvolution. It contains a_lhs, a_rhs and b_lhs. The right hand side of b we consider to be "unwrapped" and thus at
b_rhs+2pi. This yields o_lhs=max(a_lhs, b_lhs) and o_rhs=a_rhs. Now we consider the second circumvolution. It contains not only the rhs of b at b_rhs+2pi, but also the periodic repetition of a at [a_lhs+2pi, a_rhs+2pi]. This result in o_lhs=max(a_lhs+2pi, b_lhs) and o_rhs=min(a_rhs+2pi, b_rhs+2pi). Modulo 2pi shifts both down to o_lhs=a_lhs and o_rhs=min(a_rhs, b_rhs).
- case 1: both intervals spawn circumvolution one and two. The first intersection is within
[0, 4pi) the second requires the periodic repetition of one of the intervals, and thus lies within [2pi,6pi).
OLD ANSWER, deprecated:
Here my version, using numpy vector operations.
It can probably be improved by utilizing the more abstract numpy function like np.where and the like.
An other idea would be to ignore numpy and use some kind of iterator/generator function. Perhaps I will try something like that next.
import numpy as np
sectors = np.array( [[5.23,0.50], [0.7,1.8], [1.9,3.71],[4.1,5.11]] )
interval = np.array([5.7,2.15])
def normalize_sectors(sectors):
# normalize might not be the best word here
idx = sectors[...,0] > sectors[...,1]
sectors[idx,1] += 2*np.pi
return sectors
def overlapping_sectors(sectors, interval):
# 'reverse' modulo 2*pi, so that rhs is always larger than lhs"
sectors = normalize_sectors(sectors)
interval = normalize_sectors(interval.reshape(1,2)).squeeze()
# when comparing two intervals A and B, the intersection is
# [max(A.left, B.left), min(A.right, B.right)
left = np.maximum(sectors[:,0], interval[0])
right = np.minimum(sectors[:,1], interval[1])
# construct overlapping intervals
res = np.hstack([left,right]).reshape((2,-1)).T
# neither empty (lhs=rhs) nor 'reversed' lhs>rhs intervals are allowed
res = res[res[:,0] < res[:,1]]
#reapply modulo
res = res % (2*np.pi)
return res
print(overlapping_sectors(sectors, interval))
EDIT:
Here the iterator-based version. It works just as well, but appears to be somewhat numerically inferior.
def overlapping_sectors2(sectors, interval):
i_lhs, i_rhs = interval
if i_lhs>i_rhs:
i_rhs += 2*np.pi
for s_lhs, s_rhs in sectors:
if s_lhs>s_rhs:
s_rhs += 2*np.pi
o_lhs = max(s_lhs, i_lhs)
o_rhs = min(s_rhs, i_rhs)
if o_lhs < o_rhs:
yield o_lhs % (2*np.pi), o_rhs % (2*np.pi)
print(list(overlapping_sectors2(sectors, interval)))
EDIT2:
Now supports intervals that overlap in two places.
sectors = np.array( [[30,330]] )/180*np.pi
interval = np.array( [300,60] )/180*np.pi
def normalize_sectors(sectors):
# normalize might not be the best word here
idx = sectors[...,0] > sectors[...,1]
sectors[idx,1] += 2*np.pi
return sectors
def overlapping_sectors(sectors, interval):
# 'reverse' modulo 2*pi, so that rhs is always larger than lhs"
sectors = normalize_sectors(sectors)
# if interval rhs is smaller than lhs, the interval crosses 360 degrees
# and we have to consider it as two intervals
if interval[0] > interval[1]:
interval_1 = np.array([interval[0], 2*np.pi])
interval_2 = np.array([0, interval[1]])
res_1 = _overlapping_sectors(sectors, interval_1)
res_2 = _overlapping_sectors(sectors, interval_2)
res = np.vstack((res_1, res_2))
else:
res = _overlapping_sectors(sectors, interval)
#reapply modulo
res = res % (2*np.pi)
return res
def _overlapping_sectors(sector, interval):
# when comparing two intervals A and B, the intersection is
# [max(A.left, B.left), min(A.right, B.right)
left = np.maximum(sectors[:,0], interval[0])
right = np.minimum(sectors[:,1], interval[1])
# construct overlapping intervals
res = np.hstack([left,right]).reshape((2,-1)).T
# neither empty (lhs=rhs) nor 'reversed' lhs>rhs intervals are allowed
res = res[res[:,0] < res[:,1]]
return res
print(overlapping_sectors(sectors, interval)*180/np.pi)
def overlapping_sectors2(sectors, interval):
i_lhs, i_rhs = interval
for s_lhs, s_rhs in sectors:
if s_lhs>s_rhs:
s_rhs += 2*np.pi
if i_lhs > i_rhs:
o_lhs = max(s_lhs, i_lhs)
o_rhs = min(s_rhs, 2*np.pi)
if o_lhs < o_rhs:
yield o_lhs % (2*np.pi), o_rhs % (2*np.pi)
o_lhs = max(s_lhs, 0)
o_rhs = min(s_rhs, i_rhs)
if o_lhs < o_rhs:
yield o_lhs % (2*np.pi), o_rhs % (2*np.pi)
else:
o_lhs = max(s_lhs, i_lhs)
o_rhs = min(s_rhs, i_rhs)
if o_lhs < o_rhs:
yield o_lhs % (2*np.pi), o_rhs % (2*np.pi)
print(np.array(list(overlapping_sectors2(sectors, interval)))*180/np.pi)
