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.
What would it take to do this?
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.
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
As the shape is an exact rectangle, the easiest is to rotate all points by the angle
-arctan((4859302-4856938)/(677241-670388))
Doing so, the rectangle becomes axis-aligned and you just have to perform four coordinate comparisons. Rotations are easy to compute with complex numbers.
In fact you can simply represent all points as complex numbers, compute the vector defined by some side, and multiply everything by the conjugate.
A slightly different approach is to consider the change of coordinate frame that brings some corner to the origin and two incident sides to (1,0) and (0,1). This is an affine transformation. Then your test boils down to checking insideness to the unit square.