Determine if points are within a rotated rectangle (standard Python 2.7 library only) [duplicate]

Question

I have a rotated rectangle with these coordinates as vertices:

   1   670273   4879507
   2   677241   4859302
   3   670388   4856938
   4   663420   4877144

---

---

And I have points with these coordinates:

670831  4867989
675097  4869543

Using only the Python 2.7 standard library, I want to determine if the points fall within the rotated rectangle.

• I am not able to add additional Python libraries to my Jython implementation

---

What would it take to do this?

Answer 1

A line equation of the form ax+by+c==0 can be constructed from 2 points. For a given point to be inside a convex shape, we need testing whether it lies on the same side of every line defined by the shape's edges.

In pure Python code, taking care of writing the equations avoiding divisions, this could be as follows:

def is_on_right_side(x, y, xy0, xy1):
    x0, y0 = xy0
    x1, y1 = xy1
    a = float(y1 - y0)
    b = float(x0 - x1)
    c = - a*x0 - b*y0
    return a*x + b*y + c >= 0

def test_point(x, y, vertices):
    num_vert = len(vertices)
    is_right = [is_on_right_side(x, y, vertices[i], vertices[(i + 1) % num_vert]) for i in range(num_vert)]
    all_left = not any(is_right)
    all_right = all(is_right)
    return all_left or all_right

vertices = [(670273, 4879507), (677241, 4859302), (670388, 4856938), (663420, 4877144)]

The following plot tests the code visually for several shapes. Note that for shapes with horizontal and vertical lines usual line equations could provoke division by zero.

import matplotlib.pyplot as plt
import numpy as np

vertices1 = [(670273, 4879507), (677241, 4859302), (670388, 4856938), (663420, 4877144)]
vertices2 = [(680000, 4872000), (680000, 4879000), (690000, 4879000), (690000, 4872000)]
vertices3 = [(655000, 4857000), (655000, 4875000), (665000, 4857000)]
k = np.arange(6)
r = 8000
vertices4 = np.vstack([690000 + r * np.cos(k * 2 * np.pi / 6), 4863000 + r * np.sin(k * 2 * np.pi / 6)]).T
all_shapes = [vertices1, vertices2, vertices3, vertices4]

for vertices in all_shapes:
    plt.plot([x for x, y in vertices] + [vertices[0][0]], [y for x, y in vertices] + [vertices[0][1]], 'g-', lw=3)
for x, y in zip(np.random.randint(650000, 700000, 1000), np.random.randint(4855000, 4880000, 1000)):
    color = 'turquoise'
    for vertices in all_shapes:
        if test_point(x, y, vertices):
            color = 'tomato'
    plt.plot(x, y, '.', color=color)
plt.gca().set_aspect('equal')
plt.show()

PS: In case you are running a 32-bit version of numpy, with this size of integers it might be necessary to convert the values to float to avoid overflow.

If this calculation needs to happen very often, the a,b,c values can be precalculated and stored. If the direction of the edges is known, only one of all_left or all_right is needed.

When the shape is fixed, a text version of the function can be generated:

def generate_test_function(vertices, is_clockwise=True, function_name='test_function'):
    ext_vert = list(vertices) + [vertices[0]]
    unequality_sign = '>=' if is_clockwise else '<='
    print(f'def {function_name}(x, y):')
    parts = []
    for (x0, y0), (x1, y1) in zip(ext_vert[:-1], ext_vert[1:]):
        a = float(y1 - y0)
        b = float(x0 - x1)
        c = a * x0 + b * y0
        parts.append(f'({a}*x + {b}*y {unequality_sign} {c})')
    print('    return', ' and '.join(parts))

vertices = [(670273, 4879507), (677241, 4859302), (670388, 4856938), (663420, 4877144)]
generate_test_function(vertices)

This would generate a function as:

def test_function(x, y):
    return (-20205.0*x + -6968.0*y >= -47543270741.0) and (-2364.0*x + 6853.0*y >= 31699798882.0) and (20206.0*x + 6968.0*y >= 47389003912.0) and (2363.0*x + -6853.0*y >= -31855406372.0)

This function then can be copy-pasted and optimized by the Jython compiler. Note that the shape doesn't need to be rectangular. Any convex shape will do, allowing to use a tighter box.

Answer 2

Take three consequent vertices A, B, C (your 1,2,3)

Find lengths of sides AB and BC

lAB = sqrt((B.x - A.x)^2+(B.y - A.y)^2)

Get unit (normalized) direction vectors

uAB = ((B.x - A.x) / lAB, (B.y - A.y) / lAB) 

For tested point P get vector BP

BP = ((P.x - B.x), (P.y - B.y)) 

And calculate signed distances from sides to point using cross product

SignedDistABP = Cross(BP, uAB) = BP.x * uAB.y - BP.y * uAB.x
SignedDistBCP = - Cross(BP, uBC) = - BP.x * uBC.y + BP.y * uBC.x

For points inside rectangle both distances should have the same sign - either negative or positive depending on vertices order (CW or CCW), and their absolute values should not be larger than lBC and lAB correspondingly

Abs(SignedDistABP) <= lBC
Abs(SignedDistBCP) <= lAB