Fastest way to check if any point from a list is in bwteen two fixed points A and B

Question

I have two large sets of points L1 and L2 and a fixed point B. I need an efficient way to filter L1 and remove all points which their path (straight line) to B is blocked by any point from L2.

I tried the following:

# Measure distance between two points
def distance (p1, p2):
    return ((p2[1] - p1[1])**2 + (p2[0] - p1[0])**2)**0.5

# Check if point c is on same line between points a and b
def inBetween (a, b, c):
    return distance (a, b) == distance (a, c) + distance (b, c)

def isPathBlocked (L1, L2, B):
    for A in L1:
        for C in L2:
            if inBetween (A, B, C):
                print "path blocked"

Although it works but it's very CPU intensive considering that L1 and L2 can have more than 50,000 points in each and therefor the nested loops would execute 50,000 * 50,000 times. Is there a better way to do this or any good algorithm that I can use here?

Note:

I should have mentioned earlier but all points coords in L1 and L2 are integers.

Note 2:

B has integers coords and all points in L1 and L2 are within a certain radius from B (radius is about 10 to 10000).

Solutions:

Thanks to the great answers below I'm here adding an example of how to use the main two methods provided (hash and compare, sort and merge compare) with thorough explanation for those who are new to algorithms which might find it useful:

l1 = [(x,y) for x in range (10) for y in range(10)]
l2 = [(5,3), (5,1), (8,2), (4,3), (1,5), (3,5), (7,7), (3,7)]
b = (5,5)

# instead of using atan2 to calculate polar angles, to improve performance we only 
# calculate (y/x, x>0) which has a unique value for different angles
# if x == 0 (angle is eigher 90 or -90 degrees) then we use (None, y>0) instead
def toPolar (p1, p2):
    if p2[0] != p1[0]:
        return (float (p2[1] - p1[1]) / float (p2[0] - p1[0]), p2[0] < p1[0])
    if p2[1] != p1[1]:
        return (None, p2[1] < p1[1])

# Check if b is between a and c, we don't need to check if they
# are cllinear because they are already when they have same polar angle
def inBetween (a, b, c):
    if a[0] <= b[0] <= c[0] or c[0] <= b[0] <= a[0]:
        return a[1] <= b[1] <= c[1] or c[1] <= b[1] <= a[1]
    return False

# b stationary point, l1 the points of interest, l2 blockers
def hashFilter (b, l1, l2):
    hashTable = {}
    # create a hash table of l2 points using polar angles as keys,
    # if two points have the same polar angle keep only the one that
    # is closer to b
    for blocker in l2:
        key = toPolar (b, blocker)
        prev = hashTable.get (key)
        if prev == None or inBetween (b, blocker, prev):
            hashTable[key] = blocker
    unBlocked = []
    blocked = []
    # Remove the points from l1 if they are blocked by any point from l2
    for point in l1:
        key = toPolar (b, point)
        blocker = hashTable.get (key)
        if blocker == None or not inBetween (b, blocker, point):
            unBlocked.append (point)
        else: blocked.append (point)
    return sorted (blocked)

def sortFilter (b, l1, l2):
    # Convert cartesian to polar (sort of)
    points = [(toPolar(b, x), x) for x in l1]
    # Sort using polar angle as key
    points = sorted (points, key = lambda x: x[0])
    # Convert cartesian to polar (sort of)
    blockers = [(toPolar(b, x), x) for x in l2]
    # Sort using polar angle as key
    blockers = sorted (blockers, key = lambda x: x[0])
    unBlocked = []
    blocked = []
    ind = 0
    length = len (blockers)
    for point in points:
        if ind >= length: break
        isBlocked = False
        subInd = None
        # Increase index of blockers until we reach a matching polar angle
        while ind < length and blockers[ind][0] <= point[0]:
            if blockers[ind][0] == point[0]:
                # we need to check every point in points(l1) for being blocked
                # by comparing it to all points in blockers(l2) with
                # identical polar angle.
                # However because there could be multiple points and multiple
                # blockers with the same polar angle, we store the ind of the
                # first identical polar angle that appears in the blockers list,
                # so that the next item in points list can loop again over
                # this group of blockers
                if subInd == None: subInd = ind
                if inBetween (b, blockers[ind][1], point[1]):
                    blocked.append (point[1])
                    isBlocked = True
                    # If we found out that the point is blocked then we break
                    # out of the loop to test the next item in points
                    break
            ind += 1
        if not isBlocked: unBlocked.append (point[1])
        # Move the index back to the first item in the blockers group with
        # same polar angle to the last point checked
        if subInd != None: ind = subInd
    return sorted (blocked)

print hashFilter (b, l1, l2)
print sortFilter (b, l1, l2)

Answer 1

Given the constraints (small integer coordinates), we can just hash slopes instead of sorting. This depends on the facts that (1) float division is correctly rounded (2) the minimum gap between distinct fractions with denominators bounded by 10000 is at least approximately 1/10000^2, which is safely more than 1 ulp.

def key(b, p):
    if p[0] != b[0]:
        return (float(p[1] - b[1]) / float(p[0] - b[0]), p[0] < b[0])
    if p[1] != b[1]:
        return (None, p[1] < b[1])
    return None


def between1(b, p, q):
    return b <= p <= q or q <= p <= b


def between2_assume_collinear(b, p, q):
    return between1(b[0], p[0], q[0]) and between1(b[1], p[1], q[1])


def unblocked(b, l1, l2):
    inner = {}
    for p in l2:
        k = key(b, p)
        if k is None:
            return
        q = inner.get(k)
        if q is None or between2_assume_collinear(b, p, q):
            inner[k] = p
    for q in l1:
        p = inner.get(key(b, q))
        if p is None or not between2_assume_collinear(b, p, q):
            yield q

Here's some test code.

def demo(b, l2):
    n = 10
    l1 = [(x, y) for x in range(n) for y in range(n)]
    grid = [[' ' for x in range(n)] for y in range(n)]
    for q in unblocked(b, l1, l2):
        grid[q[1]][q[0]] = '*'
    grid[b[1]][b[0]] = '@'
    for p in l2:
        grid[p[1]][p[0]] = '2'
    for row in grid:
        print(''.join(row))


demo((5, 5), [(5, 3), (5, 1), (8, 2), (4, 3), (1, 5), (3, 5)])

Alternative key function that is slower but works for all integers:

import fractions


def key(b, p):
    d0 = p[0] - b[0]
    d1 = p[1] - b[1]
    g = fractions.gcd(d0, d1)
    if g == 0:
        return None
    g = abs(g)
    return (d0 // g, d1 // g)

Answer 2

this is a mini-implementation of what i mean (it is the same as idea #2 in MBo's answer):

import random
import math
from collections import namedtuple


Cartesian = namedtuple('Cartesian', ('x', 'y', 'sample'))
Polar = namedtuple('Polar', ('r', 'phi', 'sample'))


def to_polar(cartesian, origin):

    x_diff = cartesian.x - origin.x
    y_diff = cartesian.y - origin.y
    r = math.hypot(x_diff, y_diff)
    phi = math.atan2(y_diff, x_diff)

    return Polar(r=r, phi=phi, sample=cartesian.sample)


def to_cartesian(polar, origin):

    x = polar.r*math.cos(polar.phi) + origin.x
    y = polar.r*math.sin(polar.phi) + origin.y

    return Cartesian(x=x, y=y, sample=polar.sample)


random.seed(45432)  # get reproducibe results

SCALE = 100
N_SAMPLES = 5

B = Cartesian(x=SCALE*random.random(), y=SCALE*random.random(), sample='B')

L1_cartesian = [Cartesian(x=SCALE*random.random(),
                          y=SCALE*random.random(),
                          sample='L1')
                for _ in range(N_SAMPLES)]
L1_polar = [to_polar(cartesian=l1, origin=B) for l1 in L1_cartesian]

L2_cartesian = [Cartesian(x=SCALE*random.random(),
                          y=SCALE*random.random(),
                          sample='L2')
                for _ in range(N_SAMPLES)]
L2_polar = [to_polar(cartesian=l2, origin=B) for l2 in L2_cartesian]

all_polar = L1_polar + L2_polar
all_polar = sorted(all_polar, key=lambda x: x.phi)

print(all_polar)

what you get at the end is a list of all points, sorted by their angle respective to B (their origin sample is stored in all_polar.sample). you'd just have to cleverly iterate over that list; checking only angles that are 'close' (and points from different sample sets).

similar (or equal angle) then still does not mean a point is 'in between'; you'd still have to compare the distance to B (which is all_polar.r).

there is no need to join the two samples into one list; you could also iterate over one list and yield elements that are 'close' from the other list.

none of this has been throughly tested! and the comparison part is missing.

---

as you have added that your points are at integer coordinates, i suggest you do what David Eisenstat's answer suggests (this is his idea added to my code; credit should go to his answer): instead of calculating the polar coordinates, make a dictionary containing the result of this

from fractions import Fraction

def get_slope(cartesian, origin):
    x_diff = cartesian.x - origin.x
    y_diff = cartesian.y - origin.y

    slope = Fraction(y_diff, x_diff)
    x_sign = x_diff > 0

    return (slope, x_sign)

and then create a dictionary e.g. like this:

from collections import defaultdict

dct = defaultdict(list)
# add all the points in your sets (both)
dct[get_slope(cartesian, origin)].append(cartesian)

it is then trivial to compare the points for a given slope.