Search code examples
pythonalgorithmnumpygeometrysurface

How to find the closest points of a surface regarding some coordinate?

I have a bunch of points in 3d space (x,y and z) and want to find the closest point of surface existing among points. I read this solution but could not solve my issue. My surfaces are created by four points using:

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

I have a tilted (shown by black dots in my fig) and a normal (shown by yellow dots) surface among my points. Corners of the surfaces are:

surf_corners=[np.array([[1., 1., 1.], [1.5, 2., 1.], [3., 1., 3.], [3.5, 2., 3.]]),\
              np.array([[5., 1., 4.], [5., 2., 4.], [7., 1., 2.], [7., 2., 2.]])]

First array is the tilted one and second array is the normal one. These are points that I want to find the closet points of the surfaces:

surrounding_points=[np.array([[2., 1., 3.2], [3., 2., 3.]]),\
                    np.array([[6., 1., 2.5], [6., 2., 2.5], [6., 1., 3.5], [6., 2., 3.5]])]

Projection of the first array (shwon by blue squares) of surrounding_points should be found on the surface created by first arrays of the surf_corners and second one (shwon by red squares) also with the second suraface. My surface should be only among its four corners and not infinity. I mean the y value of calculated point should not be less or more than the corners. For example I prefer to maintain the y value of the surrounding_points. I mean for first point, I want the closet point of the created surface that its y value is 1. For simplicity I copied only data that their y values are equal to brders of my planes and in reality I may have other points with y values between the brders of the planes and still I want to find the projection with a fixed y value. I tried the following method to do so but it is just giving me the horizontal projection for one surface and I cannot hnadle two surfaces (when surf_corners and surrounding_points ar numpy arrays and not a list of numpy arrays):

pls=[]
for i in surf_corners:
    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(surrounding_points[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)

For more clarity, I uploaded a fig. green arrows show my desired projection point (closest point with the same y value). In advance, I do appreciate any help. enter image description here

Solution

  • I understand you need to solve the following problem:

    1. Given four 3D points forming a flat quadrilateral, determine the plane where those points lie.
    2. Project an arbitrary 3D point onto the plane.
    3. Check if the projected point lies inside the quadrilateral.

    Problem 1. Determine the plane

    I understand you are using the code from this answer, so I won't comment any further.

    Problem 2. Project a point

    You can calculate the projected point using

    projected_point = p - np.dot(p - p0, u) * u

    where

    • p is the point to project
    • p0 is a point on the plane (one of the points defining the quadrilateral, their centroid...)
    • u is the unit normal vector of the plane

    Edit:

    As long as the plane calculation method provides the plane coefficients (a,b,c,d) where (a,b,c) is the plane's normal vector, and it's also a unit vector, you can calculate the projection using:

    u = plane[:3]                      # Plane's normal unit vector
d = plane[3]                       # Signed distance plane - origin of coordinates
distance = np.dot(u, point) + d    # Signed distance point - plane
projected_point = point - distance * u

    Problem 3. Point inside quadrilateral

    I would recommend using the winding number method. And, if you know that your quadrilaterals are convex, you can take advantage of the optimizations described in the same link.