Overlapping radial intervals, python

Question

Suppose I have some intervals, in radian, which represent angular sectors on a circle. The first number is the beginning of the interval, the second number is the end of the interval. Both numbers are connected in counter-clockwise direction, i.e the trigonometric direction. I want to know their intersection with a test interval.

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])

#expected_res=function_testing_overlap_between_sectors_and_interval...
#<< expected_res = np.array[[5.7,0.50],[0.7,1.8],[1.9,2.15]]

The last sector does not overlap with the interval, so there is not intersection. For visualisation

import matplotlib
from matplotlib.patches import Wedge
import matplotlib.pyplot as plt

fig=plt.figure(figsize(10,10))
ax=fig.add_subplot(111) 
ax.set_xlim(-1.5,1.5)
ax.set_ylim(-1.5,1.5)
sectors_deg = sectors*180/pi
interval_deg = interval*180/pi

sec1 = Wedge((0,0), 1, sectors_deg[0,0], sectors_deg[0,1], color="r", alpha=0.8)
sec2 = Wedge((0,0), 1, sectors_deg[1,0], sectors_deg[1,1], color="r", alpha=0.8)
sec3 = Wedge((0,0), 1, sectors_deg[2,0], sectors_deg[2,1], color="r", alpha=0.8)
sec4 = Wedge((0,0), 1, sectors_deg[3,0], sectors_deg[3,1], color="r", alpha=0.8)
int_sec= Wedge((0,0), 1, interval_deg[0], interval_deg[1], color="b", alpha=0.2)


ax.plot(np.cos(np.linspace(0,2*pi,100)),np.sin(np.linspace(0,2*pi,100)),color='k')
ax.add_artist(sec1)
ax.add_artist(sec2)
ax.add_artist(sec3)
ax.add_artist(sec4)
ax.add_artist(int_sec)
ax.text(1.0,1.2,'red - sectors')
ax.text(1.0,1.1,'blue - interval')
ax.text(1.0,1.0,'overlap - intersect')
plt.show()

I haven't found any implemented solution to find intersection of angular sectors. Any help appreciated

Answer 1

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
  1. both wrap round
  2. only the comparison-interval wraps round
  3. only the sector-interval wraps round
  4. 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)