I have a bunch of points in 3d space (x,y and z) and want to find their perpendicular projection on a surface in python. My surface is created by four points using the following function:
PRECISION = 1e-8 # Arbitrary zero for real-world purposes
def plane_from_points(points):
centroid = np.mean(points, axis=0)
_, eigenvalues, eigenvectors = np.linalg.svd(points - centroid)
if eigenvalues[1] < PRECISION:
raise ValueError("Points are aligned, can't define a plane")
normal = eigenvectors[2]
d = -np.dot(centroid, normal)
plane = np.append(normal, d)
thickness = eigenvalues[2]
return plane, thickness
These are my input points stored as list for creating the surface:
surface_maker=[np.array([[44., 5., 25.],
[44., 25., 25.],
[59., 5., 75.],
[59., 25., 75.]])]
I want to ind the perpedicular projection of the following points in the created surface:
projection_point=[np.array([[49.9, 5., 65.],
[48., 17., 69.],
[47., 25., 71.9],
[60., 5., 39.],
[55., 15., 42.1],
[57., 25., 40.1]])]
I tried the following code to do so, but it is giving me the horizontal projection while i want to find the perpendilar projection:
pls=[]
for i in surface_maker:
i=i.tolist()
pl, thickness= plane_from_points(i)
pls.append(pl)
point_on_surf=[]
n_iter=1
for i in range (n_iter):
a, b, c, d = pls[i][0], pls[i][1], pls[i][2], pls[i][3]
def fun(row):
row[0] = -(d + b * row[1] + c * row[2]) / a # calculates new x
return row[0], row[1], row[2] # new x and old y and z
to_be_projected=[copy.copy(projection_point[i])]
new_points = np.asarray(list(map(fun, [x for point in to_be_projected for x in point])))
point_on_surf.append(new_points)
In fig I visualized what I want. I just used different colurs for points to make the figure more readable. For the upper thre point I showed arrows to visualize the projection. My code is giving me the points which are at the end of red arrows but I want to find the projection point that are perpendicular to the surface. In fact, my code is only calculating a new x for the projection_point. In the fig green arrows show the direction I want. I want to do so for all the points existing in projection_point.
In advance, I do appreciate any help.
Advertisement
Answer
One way to define a plane is by three points P, Q and R. (Four points do not necesarrily lie in the same plane, but yout four points do.) Altenatively, you can define a plane by one point P in the plane and a normal vector n, which you can determine via the cross product.
n = (Q − P) × (R – P)
Normalize the norml vector, so that you have a unit vector u of length 1:
u = n ∕ | n |
You can get the distance d of a point S to the plane from the dot product of the unit normal u with the difference vector from the point in the plane P and S:
d = (S − P) · u
The distance is signed: It is positive when S is on the side of the plane where u faces and negative when it is on the other side. (It is zero, it S is in the plane, of course.)
You can get the point S′, which is S pojected to the plane, by subtracting d · u from S:
S′ = S − d · u = S − ((S − P) · u) · u
So, now let’s put that into Python. First, Point and Vector classes. (Aren’t they the sameß You can see it that way, but I find distingishig between them useful: The difference of two Points is a Vector; a Point plus a Vector is a Point; Dot and cross products make sense for Vectors, but nor for Points. In any case I prefer to have a class with x, y and z members over tuples for spatial Vectors.)
Anyway, here goes:
class Point:
def __init__(self, x, y, z):
self.x = x
self.y = y
self.z = z
def __repr__(self):
return "Point(%g, %g, %g)" % (self.x, self.y, self.z)
def __sub__(self, other):
"""P - Q"""
if isinstance(other, Vector):
return Point(self.x - other.x,
self.y - other.y,
self.z - other.z)
return Vector(self.x - other.x,
self.y - other.y,
self.z - other.z)
class Vector:
def __init__(self, x, y, z):
self.x = x
self.y = y
self.z = z
def __repr__(self):
return "Vector(%g, %g, %g)" % (self.x, self.y, self.z)
def norm(self):
"""u / |u|"""
d = math.sqrt(self.x**2 + self.y**2 + self.z**2)
return Vector(self.x / d, self.y / d, self.z / d)
def __mul__(self, other):
"""dot product u · v or scaling x · u"""
if isinstance(other, Vector):
return (self.x * other.x
+ self.y * other.y
+ self.z * other.z)
return Vector(self.x * other,
self.y * other,
self.z * other)
def cross(a, b):
return Vector(a.y * b.z - a.z * b.y,
a.z * b.x - a.x * b.z,
a.x * b.y - a.y * b.x)
These are just the operations we need for your case. Multiplication does double duty: Multiplying two Vectors yields a scalar dot product; multiplying a Vector and a scalar number yields a scaled Vector.
Your plane, reduced to three points:
plane = (
Point(44.0, 5.0, 25.0),
Point(44.0, 25.0, 25.0),
Point(59.0, 5.0, 75.0)
)
The points you want to project:
points = [
Point(49.9, 5.0, 65.0),
Point(48.0, 17.0, 69.0),
Point(47.0, 25.0, 71.9),
Point(60.0, 5.0, 39.0),
Point(55.0, 15.0, 42.1),
Point(57.0, 25.0, 40.1)
]
And the projection code:
x = plane[1] - plane[0]
y = plane[2] - plane[0]
u = cross(x, y).norm()
for p in points:
d = (p - plane[0]) * u
pp = p - u * d
print(pp)
This will yield the points projected into the plane:
Point(55.4963, 5, 63.3211) Point(56.4404, 17, 66.4679) Point(57.156, 25, 68.8532) Point(49.1743, 5, 42.2477) Point(49.6147, 15, 43.7156) Point(49.2294, 25, 42.4312)
That plane is infinitely large, so that the projected points don’t necessarily lie between the four points of your surface.