How to split some points in 3D based on some planes in python

In 3D space I have some points (x, y and z) and want to split them by reconstructing some exisiting planes. These are my simplified points:

points=[[np.array([[20., 20., 60.], [20., 30., 65.], [55., 30., 80.], [80., 10., 60.]]),
         np.array([[20., 10., 55.], [60., 30., 70.], [70., 15., 20.]])]]

The big list has one sublist which has two arrays. I reality points has more subliststs and each sublist may also have several arrays. Then, In this simlified example I have two planes among the points. I have four corners of the two surfaces:

four_corners= [[np.array([[50., 5., 5.],
                          [50., 45., 5.],
                          [70., 45., 95.],
                          [70., 5., 95.]]),
                np.array([[30., 5., 95.],
                          [30., 45., 95.],
                          [60., 5., 5.],
                          [60., 45., 5.]])]]

four_corners has one sublist and this sublist has also two array, i.e. two surfaces. Then, each array of the first sublist of points should be devided into three parts because two surfaces are passing through them. In case of having one plane, I will have two splits for each array. I want to split my point based on the surfaces. I tried to modify this solution but I could not reach what I want. I tried the following code but it was not successful:

splitted_arr=[]
for m in points:
    for surfaces in four_corners:
        for i, j in zip (m, surfaces):
            j=np.array(j)
            corners = [j[0], j[1], j[2]]
            v1 = np.subtract(corners[1], corners[0])
            v1 = v1 / np.linalg.norm(v1)
            v2 = np.subtract(corners[-1], corners[0])
            v2 = v2 / np.linalg.norm(v2)
            orth = np.cross(v1, v2)
            for p in i:
                p_moved = np.subtract(p, corners[0])
                d = np.dot(orth, p_moved)
                if d < 0:
                    splitted_arr.append (p)
                else:
                    splitted_arr.append (p)

I want to get:

[[np.array([[20., 20., 60.], [20., 30., 65.]]), np.array([[55., 30., 80.]]), np.array([[80., 10., 60.]]),
  np.array([[20., 10., 55.]]), np.array([[60., 30., 70.]]), np.array([[70., 15., 20.]])]]

The key rule is that I want to split my points always using existing planes. In advance, I do appreciate any help and contribution.

Answer

The general principle of your approach appears sound. The fact that you call splitted_arr.append (p) in each branch of your case distinction makes that pretty much moot. So is say there is your main problem.

Some observations:

• Note that you never used j[3]. So you might simplify things by defining your planes using only three points instead of four, and in doing so eliminate the possibility of inconsistency where the for points don’t lie in a single plane.

• You might want to compute ortho for each surface outside the loop over points, to reduce repeated computations and improve performance.

• The normalisation step where you divide each of v1, v2 by its length is not needed for your application so you can omit that. This will impact the magnitude of the result but not the sign.

• Instead of subtracting corners[0] before the for product, you may as well subtract the for product of that and p (i.e. np.dot(orth, corners[0])) from the result of the for product. This replaces a vector subtraction by a scalar subtraction. Or, even better, instead of subtracting that number and comparing against zero, compare against that number right away. Again that dot product of the corner is something you can compute outside the loop over all the points. All of this will help performance without changing the result (apart from potential rounding differences).

Precomputing the ortho vector and corner dot product for each plane is conceptually the same as converting it to normal form, i.e. finding a, b, c, d such that a*x+b*y+c*z=d for each point on the plane. That formulation might be particularly attractive for readers who are not building on the numpy library.

In your question you write that (my reformatting aside) you want to get

[[np.array([[20., 20., 60.], [20., 30., 65.]]),
  np.array([[55., 30., 80.]]),
  np.array([[80., 10., 60.]]),
  np.array([[20., 10., 55.]]),
  np.array([[60., 30., 70.]]),
  np.array([[70., 15., 20.]])]]

Why? Why are the first two points in a single group, and all the others separate? And do you really want to combine multiple points into a single 2D array again?

In a comment below you write that expressing the relationship relative to several planes is part of what makes this difficult. In essence given a point and a list of (oriented) planes, you can decide whether the point is on the positive or on the negative side of each of these planes (ignoring the exactly-on-the-plane case which for rounding reasons shouldn’t be relied on anyway). So for each point and plane you get one bit of relative information, and for a list of such planes you get a list of these comparison results. You can use them as a key to group points that are not separated from one another by any of the planes. The choice which side of each plane is positive and which side is negative doesn’t matter too much, as long as it’s done consistently i.e. each ortho gets computed only once or at least always in the same way.

https://ideone.com/oo5GFg has a working example putting it all together.