how to pass a surface through four points and find the projection of adjacent points on the surface

Question

I have some kind of a complicated geometrical question. I have two series of points represented by their x,y and z as numpy arrays (in reality I have thousands in each series but make it here simple for a better understanding). First data set (arr_1) represents a cutting surface. This surface is usually vertical and displaces horizontal surfaces. The second set is points of the vertical surfaces that are displaces by a surface (the surface that arr_1 is representing). At the moment I have some redundant data in these two sets. Firstly, in arr_1 I prefer to have only four corners points and make a surface by myself using them. Secondly, in arr_2 I also do not need the data that are too close to arr_1. Finally, I want to terminate the points excising in both sides of the surface which is made by four corners of the arr_1.

import numpy as np
arr_1= np.array ([[31.95548952,  7.5       , 12.5       ],
       [31.95548952, 22.5       , 12.5       ],
       [37.5       ,  7.5       , 24.20636043],
       [37.5       , 22.5       , 24.20636043],
       [43.78154278,  7.5       , 37.5       ],
       [43.78154278, 22.5       , 37.5       ],
       [55.32209575,  7.5       , 62.5       ],
       [55.32209575, 22.5       , 62.5       ],
       [62.5       ,  7.5       , 78.50696445],
       [62.5       , 22.5       , 78.50696445],
       [66.52446985,  7.5       , 87.5       ],
       [66.52446985, 22.5       , 87.5       ]])
arr_2= np.array([[87.5       ,  7.5       , 49.99997914],
       [87.5       , 22.5       , 50.00001192],
       [62.5       ,  7.5       , 50.00004172],
       [62.5       , 22.5       , 50.00007749],
       [46.8747884 ,  7.5       , 62.5       ],
       [46.87483609, 22.5       , 62.5       ],
       [37.5       ,  7.5       , 69.99973059],
       [37.5       , 22.5       , 69.99977231],
       [12.5       ,  7.5       , 70.00012398],
       [12.5       , 22.5       , 70.00015974]])

Then, I find the distances between each set of these two arrays using the following code:

from scipy.spatial import distance
distances=distance.cdist(arr_1, arr_2)

When I run it, distances is an array with 12 rows and 10 columns. Now, I want to remove the points of arr_2 which closer than a threshold (let's say 10) to any point of arr_1. I showed these two points in a black rectangle in the first uploaded photo. A good solution is proposed here, but it does not solve my problem because I want to compare the distance of each row of arr_1 with each row of arr_2. I appreciate any solution to do it. In fact, arr_1 contains the points of a surface but I do not need all of them to make my surface. I only pick for corners of this set to make my surface. I use a highly time consuming for loop, so I appreciate any faster method to find corners of my points:

corner1 = np.array(arr_1.min (axis=0)) # this gives the row containing MIN values of x,y and z
corner2 = np.array([])
corner4 = np.array([])
corner3 = np.array(arr_1.max (axis=0)) # this gives the row containing MAX values of x,y and z
# the next block will find the corner in which x and y are minimum and z is maximum
for i in arr_1[:,0]:
    if i == max (arr_1[:,0]):
        for j in arr_1[:,1]: 
            if j == min (arr_1[:,1]):
                for h in arr_1[:,2]:
                    if h == max (arr_1[:,2]):
                        corner2 = np.append(corner2, np.array([i,j,h]))
                        corner2=corner2[0:3]
# the next block will find the corner in which x and z are minimum and y is maximum
for m in arr_1[:,0]:
    if m == min (arr_1[:,0]):
        for n in arr_1[:,1]: 
            if n == max (arr_1[:,1]):
                for o in arr_1[:,2]:
                    if o == min (arr_1[:,2]):
                        corner4 = np.append(corner4, np.array([m,n,o]))
                        corner4=corner4[0:3]

Finally, after extracting four corners, I want to make a surface using them, and find the perpendicular (green arrows) or horizonal (red arrows) projection of adjacent points of arr_2 on it. I have no idea idea how to find the surface and projections. Thanks for reading and payying attention to my detailed issue! I appreciate if anyone propose any solution to any part of it.

Answer 1

Let's split this problem up into a couple of parts.

1. You have a bunch of datapoints that describe a noisy plane, arr_1.
2. You want to know where it intersects another plane, described by arr_2.
3. You want to establish some threshold in arr_2 about that intersection.

The approach I show here is to assume that the data is a measurement of some truth value, and you want to perform these operations on your best guess at that value, not the raw data. To that end:

Part 1: Least squares fit to plane

There are a couple of different ways to describe a plane, such as a normal vector and a point. The simplest for least-squares fitting is probably

a * x + b * y + c * z = 1

Assuming that your data is reasonably well behaved, there should be no catch to doing a simple fit using the equation

arr_1 @ [[a], [b], [c]] = 1   # almost python pseudo code

Since there is no single solution possible to this with more than four points, you can run np.linalg.lstsq to get the values optimized according to the MSE:

plane_1 = np.linalg.lstsq(arr_1, np.ones((arr_1.shape[0], 1), dtype=arr_1.dtype))[0].ravel()

If arr_2 is a plane as well, you do the same thing to it to get plane_2.

Part 2: Intersection of planes

Many solutions for the intersection of planes rely on the normal vectors of the planes. I will assume that both planes have dependence on the Y coordinate (which seems safe from the plots). In that case, you can follow this answer on Math Stack Exchange. Setting y = t, you can solve for the line from the system

a1 * x + c1 * z = 1 - b1 * t
a2 * x + c2 * z = 1 - b2 * t

Here, the vector [a1, b1, c1] is plane_1. After working out the details, you get

m = np.cross(plane_1, plane_2)
b = np.array([plane_1[2] - plane_2[2], 0, plane_2[0] - plane_1[0]]) / m[1]
m /= m[1]
line = (m * t + b)

This is a parametrization for any value of t.

Part 3: Distance from point to line

To threshold the values of arr_2 about the line computed above in terms of m and b, you need a formula for the distance between a point and a line. This can be done with e.g., the method in the posts here and here.

For example, a single point p can be processed as follows:

t = (p - b).dot(m) / m.dot(m)
dist = np.sqrt(np.sum((p - (m * t + b))**2))

If you are only interested in doing the thresholding, you can compare dist**2 to the square of the threshold and save some cycles on the square root since both functions are monotonic.

TL;DR

Inputs: arr_1, arr_2, distance_threshold.

# Find planes in data
plane_1 = np.linalg.lstsq(arr_1, np.ones((arr_1.shape[0], 1), dtype=arr_1.dtype))[0].ravel()
plane_2 = np.linalg.lstsq(arr_2, np.ones((arr_2.shape[0], 1), dtype=arr_2.dtype))[0].ravel()

# Find intersection line, assume plane is not y=const
m = np.cross(plane_1, plane_2)
b = np.array([plane_1[2] - plane_2[2], 0, plane_2[0] - plane_1[0]]) / m[1]
m /= m[1]

# Find mask for arr_2
t = ((arr_2 - b).dot(m) / m.dot(m))[:, None]
dist2 = np.sum((arr_2 - (m * t + b))**2, axis=1)
mask = dist2 >= distance_threshold**2

# Apply mask
subset = arr_2[mask, :]

Appendix 1: RMSE

The really nice thing about this approach (besides the fact that it limits your algorithm to ~O(n)), as mentioned before, is that it de-noises your data. You can use the result of the least squares fits to compute the RMSE of the plane data about the fit to get a sense of just how planar your measurements really are. The second return value of lstsq is the RMSE metric.

Appendix 2: Distance of point to plane

If the data in arr_2 is indeed not planar, you can subset it a little differently. Rather than finding the intersection of a pair of planes, you can use the formula for the distance between a single point and a plane directly, as given here:

np.abs(p * plane_1 - 1) / np.sqrt(plane1.dot(plane_1))

The code then becomes

# Find planes in data
plane_1 = np.linalg.lstsq(arr_1, np.ones((arr_1.shape[0], 1), dtype=arr_1.dtype))[0].ravel()
plane_2 = np.linalg.lstsq(arr_2, np.ones((arr_2.shape[0], 1), dtype=arr_2.dtype))[0].ravel()

# Find mask for arr_2
dist2 = (arr_2 * plane_1 - 1)**2 / plane_1.dot(plane_1)
mask = dist2 >= distance_threshold**2

# Apply mask
subset = arr_2[mask, :]

Answer 2

Here's some code I put together - it covers the first part of your question - which is to find and define "surfaces" from your collections of points.

Any flat surface (a plane) can be defined by the formula ax + by + cz + d = 0 where a,b,c are coefficients and d is some constant. We can write that into a function like the one below:

def linear_plane(data, a, b, c):
    x = data[0]
    y = data[1]
    z = data[2]
    return (a * x) + (b * y) + (c*z)

Then, after some imports, we can use scipy's curve_fit to search for and find the parameters that best fit that function for each of your arrays. curve_fit takes input values does the linear_plane function and tries to tweak a,b,c so that the results of that function match the a list containing lots of -1's. This comes from rearranging the equation for the case where d = 1

from scipy.optimize import curve_fit

v1,err1 = curve_fit(linear_plane,arr_1.T,np.repeat(-1,len(arr_1)))
v2,err2 = curve_fit(linear_plane,arr_2.T,np.repeat(-1,len(arr_2)))

Where each of v1, v2 are the coefficients a,b,c that most closely define the two planes you're talking about. err1 and err2 show the residual "errors" so watch out if these get too large.

Now the coefficients are found, you can use them to visualise the two planes:

%matplotlib inline
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter3D(xs=arr_1[:,0], ys=arr_1[:,1], zs=arr_1[:,2], c="blue")
ax.scatter3D(xs=arr_2[:,0], ys=arr_2[:,1], zs=arr_2[:,2], c="red")
xx1, yy1 = np.meshgrid([min(arr_1[:,0]),max(arr_1[:,0])],[min(arr_1[:,1]),max(arr_1[:,1])])
xx2, yy2 = np.meshgrid([min(arr_2[:,0]),max(arr_2[:,0])], [min(arr_2[:,1]),max(arr_2[:,1])])

# Use the coefficients in v1[0], v1[1], v1[2] to calculate z-values for all xx/yy values
z1 = (-v1[0] * xx1 - v1[1] * yy1 - 1) * 1. /v1[2]


# Use the coefficients in v2[0], v2[1], v2[2] to calculate z-values for all xx/yy values
z2 = (-v2[0] * xx2 - v2[1] * yy2 - 1) * 1. /v2[2]

# Using xx,yy and z values, plot the surfaces fit into the graph.
ax.plot_surface(xx1, yy1, z1, alpha=0.2, color=[0,0,1])
ax.plot_surface(xx2, yy2, z2, alpha=0.2, color=[1,0,0])

Here, we calculate a collection of z values for all x's and y's falling within the max and min bounds of each plane, and then plot these as surfaces showing the intersection:

I notice Mad Physicist has just posted a more complete answer, so will leave this here for now. One thing I'd like to be able to do is extend the surface fitting so that the planes don't have to be linear - like an analogous to some 3d polynomial, which should just be a case of substituting the function used in the curve_fit call.

The other part missing from here is how to calculate the line equation that describes where these two planes cross, which I believe is found in the other (better) answer.