how to divide a numpy array of coordinates using a surface

Question

I have two numpy arrays both having coordinates of some points. One array has the coordinates of a cutting planar surface (cutting_surf) and another one has coordinates of some points (points):

cutting_surf=np.array([[3., 1., 3.], [2., 1., 1.],[3., 2., 3.],\
                       [2., 2., 1.], [3., 3., 3.], [2., 3., 1.]])
points=np.array([[1.2, 3., 3.], [4., 1., 2.], [2.2, 2., 1.], [1.2, 1.5, 3.], [2.5, 1.5, 3.],\
                 [2.9, 1., 3.], [2.9, 2., 2.9], [4., 3., 1.9], [3.2, 1.2, 3.], [2.2, 3., 3.5],\
                 [2.5, 3., 3.], [2.2, 1.5, 3.5]])

Then, I want to remove some redundant coordinates from points. These coordinates are too close to the coordinates from cutting_surf. I used the following code to do so:

to_be_removed=points[np.where(np.min(distance.cdist(cutting_surf, points),axis=0)<0.5)[0],:]
cleaned_result= npi.difference(points, to_be_removed) # it removes that close points

After that, my final plan is to divid my points array into two ones. In fact, I want to really cut points array using cutting_surf. I mean I want to get my cleaned_result as:

[np.array([[2.5, 3., 3.],
          [2.5, 1.5, 3.],
          [1.2, 3., 3.],
          [1.2, 1.5, 3.],
          [2.2, 3., 3.5],
          [2.2, 1.5, 3.5]])
 np.array([[4. , 3. , 1.9],
           [4. , 1. , 2. ]])]

My fig shows the distribution of coordinates that I have. I presented just some simple coordinates here and by sorting them based on their x values I can divide the array into two one but in reality it is much more complicates. I think the only way is creating the surface using cutting_surf coordinates and then separating two clusters. I tried creating the surface using four corners of cutting_surf but I have no idea how to cluster my points using this surface:

def plane_from_points(cutting_surf):
    centroid = np.mean(cutting_surf, axis=0)
    _, eigenvalues, eigenvectors = np.linalg.svd(cutting_surf - 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

In advance, I appreciate any help and contribution.

Answer 1

If I understand correctly, you have already done the filtering by distance, so I will focus on the aspect of separating the "clean" points.

Please note: In your example, the six cutting_surf points are on two parallel lines, hence they span a plane - if your surface is not a simple plane, this simple approach will not work.

Geometrically speaking:

1. We determine a vector orth orthogonal to the cutting_surf plane.
2. For each point p of the points, we determine the dot product between orth and p relative to some point on the plane. The sign of the dot product is the "side" of the plane: All points with a positive sign are on one side, all points with a negative sign on the other. (And those with a zero dot product are on the plane.)

For the first part, we describe the cutting_surf plane in terms of two normed vectors v1 and v2, and the orthogonal vector o is just the cross product of these. This is basically an adaptation of another Stackoverflow answer:

import numpy as np
# surface points slightly reordered
cutting_surf = np.array([[3., 1., 3.], [3., 2., 3.], [3., 3., 3.],
                         [2., 1., 1.], [2., 2., 1.], [2., 3., 1.]])

# take three points, not all on the same line
# -> two vectors for the plane
corner = cutting_surf[0]
v1 = np.subtract(cutting_surf[1], corner)
v1 = v1 / np.linalg.norm(v1)

v2 = np.subtract(cutting_surf[3], corner)
v2 = v2 / np.linalg.norm(v2)

orth = np.cross(v1, v2)

For the second part, take each "clean" point, determine its difference to some point on the plane, and group by the sign of the dot product with orth:

cleaned_points = np.array([[2.5, 3., 3.], [2.5, 1.5, 3.], [1.2, 3., 3.],
      [1.2, 1.5, 3.], [2.2, 3., 3.5], [2.2, 1.5, 3.5], [4. , 3. , 1.9],
      [4. , 1. , 2. ]])

one, other = [], []
for p in cleaned_points:
    # point relative to our corner:
    p_moved = np.subtract(p, corner)
    d = np.dot(orth, p_moved)
    # in case of doubt, check for d == 0 here
    (one if d < 0 else other).append(p)
print(one)
print(other)

The output of the sides one and other seems to be the desired grouping:

[array([4. , 3. , 1.9]), array([4., 1., 2.])]
[array([2.5, 3. , 3. ]), array([2.5, 1.5, 3. ]), array([1.2, 3. , 3. ]), array([1.2, 1.5, 3. ]), array([2.2, 3. , 3.5]), array([2.2, 1.5, 3.5])]

Here is a script that renders a 3d plot of the different elements:

from matplotlib import pyplot as plt
import numpy as np

# surface points slightly reordered
cutting_surf = np.array([[3., 1., 3.], [3., 2., 3.], [3., 3., 3.],
                        [2., 1., 1.], [2., 2., 1.], [2., 3., 1.]])

# take three points, not all on the same line
# -> two vectors for the plane
corner = cutting_surf[0]
v1 = np.subtract(cutting_surf[1], corner)
v1 = v1 / np.linalg.norm(v1)

v2 = np.subtract(cutting_surf[3], corner)
v2 = v2 / np.linalg.norm(v2)

orth = np.cross(v1, v2)

fig = plt.figure()
ax = plt.axes(projection='3d')
ax.set_xlim((0.5, 3.5))
ax.set_ylim((0.5, 3.5))
ax.set_zlim((0.5, 3.5))

# plot the surface points
xs, ys, zs = ([p[i] for p in cutting_surf] for i in range(3))
ax.scatter3D(xs, ys, zs, color="grey", marker="+")

# our two plane vectors:
ax.plot([xs[0], xs[1]], [ys[0], ys[1]], [zs[0], zs[1]], color="grey", alpha=0.6)
ax.plot([xs[0], xs[3]], [ys[0], ys[3]], [zs[0], zs[3]], color="grey", alpha=0.6)
# the orthogonal vector, based on the corner point, in green
ax.plot([corner[0], corner[0] + orth[0]],
        [corner[1], corner[1] + orth[1]],
        [corner[2], corner[2] + orth[2]], color="green", alpha=0.6)

cleaned_points = np.array([[2.5, 3., 3.], [2.5, 1.5, 3.], [1.2, 3., 3.],
          [1.2, 1.5, 3.], [2.2, 3., 3.5], [2.2, 1.5, 3.5], [4. , 3. , 1.9],
          [4. , 1. , 2. ]])

for p in cleaned_points:
    # point relative to our corner:
    p_moved = np.subtract(p, corner)
    d = np.dot(orth, p_moved)
    ax.plot(*p, color="cyan" if d < 0 else "red", marker="D")
plt.show()

So here you can see the plane (grey lines), the orthogonal vector (green), and the separated "clean" points (cyan and red):

As a possible bonus, if you normalize the orthogonal vector, the dot product gives you the distance of the points to the plane, so you could merge the cleaning step into this one (see, e.g., https://mathinsight.org/distance_point_plane).