I have a 3D matrix in python:
mat = np.zeros((512,512,480))
The index of each element is the coordinates of each point in space.
And a plane:
center = np.array([251.7, 150.2, 200.1])
normal = np.array([-0.07017233, -0.99284769, 0.09658831])
I want to have a mask of mat. The value of each element is 1 if the corresponding point is on the same side as the direction of the normal, otherwise, set the value to 0.
Let n = (nx, ny, nz) be your normal vector and c = (xo,yo,zo) a point of the plane.
And given a point P whose position vector is r = (x,y,z).
If the dot product of (r-c) and n be positive, the point P is above the plane, otherwise, the point is below the plane.
You can check this by analytical geometry. But this is just the math, one way to implement this is by doing a loop through the indices of the matrix, making a vector of them, subtracting the center, and then using np.dot() to make a dot product with normal inside an if statement that will decide if the value will be 0 or 1.
Ok, this is the most efficient solution that I could think:
center = np.array([251.7, 150.2, 200.1])
normal = np.array([-0.07017233, -0.99284769, 0.09658831])
x = np.arange(512)
y = np.arange(512)
z = np.arange(480)
xmesh, ymesh, zmesh = np.meshgrid(x,y,z)
xm = xmesh - center[0]
ym = ymesh - center[1]
zm = zmesh - center[2]
p = xm*normal[0] + ym*normal[1] + zm*normal[2]
mat = np.where(p>0,1,0)
This follows the same mathematical thinking that a described in the other answer, but, with the meshgrids and the np.where command we can check the location of the point in relation to the plane in parallel.
I implemented a version with 3 four loops and it runs in more than 6 minutes, this version runs in 30 seconds. I hope that it can help you.