I'm trying to speed up a simple interpolation algorithm. The problem is as follows:
There are 4 three-dimensional points which form a rectangle in XY plane. Their Z components are arbitrary.
The 4 points form a surface between them and their midpoint, like this:
I need to project a point on that surface by the XY coordinates and getting the Z coordinate on the surface as the result:
My solution (slow):
def project_point(p1, p2, p3, test):
"""
Define a plane by 3 points and project the test point on it
"""
v1 = p3 - p1
v2 = p2 - p1
# The cross product is a vector normal to the plane
cp = np.cross(v1, v2)
a, b, c = cp
# This evaluates a * x3 + b * y3 + c * z3 which equals d
d = np.dot(cp, p3)
# Z = (d - a * X - b * Y) / c
return (d - a * test[0] - b * test[1]) / c
# Determine one of the 4 triangles and project on it
# points is a list of 3d numpy arrays
# xy is a 2d numpy array
p1, p2, p3, p4 = tuple(points)
mid_point = (p1 + p2 + p3 + p4) * 0.25
tri_index = get_triangle_index(p1, p2, p3, p4, xy) # Returns a triangle index at the coordinates xy
z = 0
if tri_index == 0:
z = project_point(p1, p3, mid_point, xy)
elif tri_index == 1:
z = project_point(p3, p4, mid_point, xy)
elif tri_index == 2:
z = project_point(p4, p2, mid_point, xy)
elif tri_index == 3:
z = project_point(p1, p2, mid_point, xy)
The provided solution works but it's slow and something tells me this has been solved before more elegantly. Can anyone point me in the right direction?
