Search code examples
pythonnumpygeometrysurface

how to find the coordinate of points projection on a planar surface

Hope doing well. I have two numpy array, both are some points in the space. Using python, I want to firstly find the surface passing the first data set (surface_maker) and then find the x,y and z of the projection adjacent opoints of the second array (contact_maker) on the created surface. surface_maker always created planar surfaces. For projection, I only want a vertical going from adjacent point toward the surface. In reality I have lots of points in both sets but I copies a simple case here:

surface_maker=np.array([[50., 15., 46.04750574],
                        [50., 5., 45.56400925],
                        [44.83018398, 5., 25.],
                        [44.76296902, 15., 25.],
                        [50., 25., 45.56400925],
                        [44.83018398, 25., 25.],
                        [59.8336792, 5., 75.],
                        [59.71483707, 15., 75.],
                        [59.8336792, 25., 75.]])
contact_maker=np.array([[10.,  5., 70.00014782],
                        [10., 15., 70.00018358],
                        [10., 25., 70.0001955 ],
                        [30.,  5., 69.99981105],
                        [30., 15., 69.99982297],
                        [30., 25., 69.99985874],
                        [70., 5., 50.00000298],
                        [70., 15., 50.00002682],
                        [70., 25., 50.00005066],
                        [90., 5., 49.99996871],
                        [90., 15., 49.99999255],
                        [90., 25., 50.00001788]])

I have tried several solutions like 1, 2 and so on. But I was successful to solve my issue. For me it is important to have the location of projection as x, y and z. The figure also shows what I want (as it shows, I need only location six adjacent point of the contact_maker projected on the surface created by surface_maker):

enter image description here

In advance, I truely appreciate any help.

Solution

  • I understand you need to solve two problems:

    • Find the plane that fits a collection of points
    • Project a second collection of points onto that plane along a specific direction

    The second problem has been fully addressed in another answer, so I'm contributing a more generic approach to the first problem.

    It's true that when you positively know that all your points lie on a plane, you may just select three non-aligned ones and calculate the plane. But your points may come from real measurements with some noise, and you may wish to find the plane that best fists your points.

    The following function solves the general problem of finding the plane that best fits a collection of points. See the explanations in the comments:

    import numpy as np
PRECISION = 1e-8    # Arbitrary zero for real-world purposes

def plane_from_points(points):
    # The adjusted plane crosses the centroid of the point collection
    centroid = np.mean(points, axis=0)

    # Use SVD to calculate the principal axes of the point collection
    # (eigenvectors) and their relative size (eigenvalues)
    _, values, vectors = np.linalg.svd(points - centroid)

    # Each singular value is paired with its vector and they are sorted from
    # largest to smallest value.
    # The adjusted plane plane must contain the eigenvectors corresponding to
    # the two largest eigenvalues. If only one eigenvector is different
    # from zero, then points are aligned and they don't define a plane.
    if values[1] < PRECISION:
        raise ValueError("Points are aligned, can't define a plane")

    # So the plane normal is the eigenvector with the smallest eigenvalue
    normal = vectors[2]

    # Calculate the coefficients (a,b,c,d) of the plane's equation ax+by+cz+d=0.
    # The first three coefficients are given by the normal, and the fourth
    # one (d) is the plane's signed distance to the origin of coordinates
    d = -np.dot(centroid, normal)
    plane = np.append(normal, d)

    # If the smallest eigenvector is close to zero, the collection of
    # points is perfectly flat. The larger the eigenvector, the less flat.
    # You may wish to know this.
    thickness = values[2]

    return plane, thickness

    You can check this:

    >>> surface_maker=np.array([[50., 15., 46.04750574], [50., 5., 45.56400925], [44.83018398, 5., 25.], [44.76296902, 15., 25.], [50., 25., 45.56400925], [44.83018398, 25., 25.], [59.8336792, 5., 75.], [59.71483707, 15., 75.], [59.8336792, 25., 75.]])
>>> plane, thickness = plane_from_points(surface_maker)
>>> print(plane)
[-0.95725318  0.          0.28925136 35.2806339 ]
>>> print(thickness)
1.3825669490602308

    So, in fact, your point distribution is not flat (thickness clearly different from zero), and you can't just select three arbitrary points to solve your problem.