Python point lookup (coordinate binning?)

Greetings,

I am trying to bin an array of points (x, y) into an array of boxes [(x0, y0), (x1, y0), (x0, y1), (x1, y1)] (tuples are the corner points)

So far I have the following routine:

def isInside(self, point, x0, x1, y0, y1):
    pr1 = getProduct(point, (x0, y0), (x1, y0))
    if pr1 >= 0:
        pr2 = getProduct(point, (x1, y0), (x1, y1))
        if pr2 >= 0:
            pr3 = getProduct(point, (x1, y1), (x0, y1))
        开发者_开发问答    if pr3 >= 0:
                pr4 = getProduct(point, (x0, y1), (x0, y0))
                if pr4 >= 0:
                    return True
    return False

def getProduct(origin, pointA, pointB):
    product = (pointA[0] - origin[0])*(pointB[1] - origin[1]) - (pointB[0] - origin[0])*(pointA[1] - origin[1])
    return product

Is there any better way then point-by-point lookup? Maybe some not-obvious numpy routine?

Thank you!

If I understand your problem correctly then the following should work assuming that your points are also 2-tuples.

def in_bin(point, lower_corner, upper_corner):
    """
    lower_corner is a 2-tuple - the coords of the lower left hand corner of the
    bin.
    upper_corner is a 2-tuple - the coords of the upper right hand corner of the
    bin.
    """
    return lower_corner <= point <= upper_corner

if __name__ == '__main__':
    p_min = (1, 1) # lower left corner of bin
    p_max = (5, 5) # upper right corner of bin

    p1 = (3, 3) # inside
    p2 = (1, 0) # outside
    p3 = (5, 6) # outside
    p4 = (1, 5) # inside

    points = [p1, p2, p3, p4]

    for p in points:
        print '%s in bin: %s' % (p, in_bin(p, x_min, x_max))

This code shows that you can compare tuples directly - there is some information in the documentation about this: http://docs.python.org/tutorial/datastructures.html#comparing-sequences-and-other-types

Are these boxes axis aligned? i.e. are the edges parallel to the coordinate axes? If so this can be done quite efficiently with vectorized comparisons on NumPy arrays.

def in_box(X, B):
    """
    Takes an Nx2 NumPy array of points and a 4x2 NumPy array of corners that 
    form an axis aligned box.
    """
    xmin = B[:,0].min(); xmax = B[:,0].max()
    ymin = X[:,1].min(); ymax = X[:,1].max()
    return X[X[:,0] > xmin & X[:,0] < xmax & X[:,1] > ymin & X[:,1] < ymax]

amending to >= and <= if you prefer them to be inclusive.

If you need it for an arbitrary quadrilateral, matplotlib actually has a routine matplotlib.nxutils.points_inside_poly that you could use (if you have it installed) or else copy it (it's BSD-licensed). See this page for a discussion of the algorithms used and other algorithms for inside-a-polygon tests.

Without too much change, your code can be compacted down to:

def isInside(self, point, x0, x1, y0, y1):
    return getProduct(point, (x0, y0), (x1, y0)) >= 0 and
           getProduct(point, (x1, y0), (x1, y1)) >= 0 and
           getProduct(point, (x1, y1), (x0, y1)) >= 0 and
           getProduct(point, (x0, y1), (x0, y0)) >= 0

def getProduct(origin, pointA, pointB):
    product = (pointA[0] - origin[0])*(pointB[1] - origin[1]) - (pointB[0] - origin[0])*(pointA[1] - origin[1])
    return product

Your solution is O(N) where N is number of points. If N is large enough and you are running the query isInside a lot of times, you might considering sorting the points and then using binary search in order to find the relevant points.

As always, first profile whether you really need this optimisation.

I used a similar routine to do colourmapped density plots:

#calculate densities
rho = zeros((nx,ny));
for i in range(N):
    x_sample = int(round(ix[i]))
    y_sample = int(round(iy[i]))

    if (x_sample > 0) and (y_sample > 0) and (x_sample<nx) and (y_sample<ny):
        rho[y_sample,x_sample] = rho[y_sample,x_sample] + 1

Instead of counting density you can store the x and y samples.

If you really do need to use getProduct... packing, unpacking and good variable names ftw!

def isInside(self, point, x0, x1, y0, y1):
    A = x0,y0
    B = x1,y0
    C = x1,y1
    D = x0,y1

    return getProduct(point, A, B) and
           getProduct(point, B, C) and
           getProduct(point, C, D) and
           getProduct(point, D, A)

def getProduct(origin, pointA, pointB):
    xA,yA = pointA
    xB,yB = pointB
    x,y = point

    return (xA - x)*(yB - y) - (xB - x)*(yB - y)

Are you sure you need such a complicated check to begin with?

def isInside(self, point, x0, y0, x1, y1):
  x,y = point
  if x0 > x1: x0,x1 = x1,x0 #these cause no
  if y0 > y1: y0,y1 = y1,y0 #side effect.

  return x0 <= x <= x1 and y0 <= y <= y1

Assuming that your boxes are rectangular, do not overlap, and have no gaps, then why don't you just call numpy.histogram2d? See the numpy docs.