3D hyper-plane in a numpy array

Question

I have the vertices of a hyper-plane like so: vertice_list = [[0, 0, 145], [0, 229, 145], [251, 0, 157], [251, 229, 157]]. I need a numpy array of shape matrix_shape = (251, 229, 388), that has a value of 1 on the hyper-plane and 0 otherwise. How can I make a function which is input a vertice_list and a matrix_shape that returns my desired matrix? (please check the image for clarification, suppose the blue hyperplane has a value of 1 and the space out of the hyperplane has a value of zero)

Your `matrix_shape` is 3D, so a hyperplane for the space it defines would be 2D? But the hyperplane you define is 3D, so how is `matrix_shape` oriented in your 4D space? Or is it mission a dimension? — Grismar, Mar 05 '22 at 21:37
OK, looks like I misread, because although you provide four coordinates (three would suffice), they do define a 2D plane in 3D space. You ask "array [..] that has a value of 1 on the hyper-plane and 0 otherwise" - do you mean 'each cell that is intersected by the hyperplane should be on 1 and cell that aren't intersected by it should be 0'? — Grismar, Mar 05 '22 at 22:07

score 0 · Answer 1 · answered Mar 06 '22 at 09:09

I think the right way is probably to extend Bresenham's algorithm to 3 dimensions; a quick google search found this paper, which I don't have access to, and which looks quite complicated.

A naive way would be to repeat Bresenham's algorithm per 2D slice of the 3D space. I think you'd have to pick the direction of iteration based on slope (normal) of the plane, a bit like one would when picking to increment x or y by 1 when doing Bresenham in 2D.

There's a much simpler, but probably less accurate method:

Find n and c for equation of the plane dot(n, x) = c from the vertices
Use a simple condition based on this to set the elements of your 3D array. This is demonstrated below for random points in a small space; I tried your example but displaying it with Matplotlib's voxels took way too long (I killed it after about a minute).

import numpy as np
import matplotlib.pyplot as plt

def get_plane_parameters(points):
    """Find normal and offset of plane fitting points"""
    # https://stackoverflow.com/a/18789789/1008142
    # I *think* that answer has a mistake in using the last column
    # of V, rather than the last row.
    cm = np.mean(points, axis=0) # column mean
    points0 = points - cm[np.newaxis]
    u, s, vh = np.linalg.svd(points0)
    # todo: use s to check collinearity, or non-coplanarity
    n = vh[-1]
    c = n @ cm
    return n, c

# limits on dimensions
# ax.voxels is slow!  Use something else for bigger cases.
xdim = 20
ydim = 20
zdim = 20

# simpler test cases
# x0, y0, z0 = [0,   0, 10]
# x1, y1, z1 = [0,  10, 10]
# x2, y2, z2 = [10, 10, 10]

# x0, y0, z0 = [10, 0,  0]
# x1, y1, z1 = [10, 0,  10]
# x2, y2, z2 = [10, 10, 10]

# x0, y0, z0 = [0, 0, 0]
# x1, y1, z1 = [0, 10, 10]
# x2, y2, z2 = [10, 10, 10]

# random points
x0, x1, x2 = np.random.choice(np.arange(xdim), size=3)
y0, y1, y2 = np.random.choice(np.arange(ydim), size=3)
z0, z1, z2 = np.random.choice(np.arange(zdim), size=3)

vertice_list = np.array([[x0, y0, z0],
                         [x1, y1, z1],
                         [x2, y2, z2]])
n, c = get_plane_parameters(vertice_list)

if not np.max(abs(vertice_list @ n - c)) < 1e-7:
    raise ValueError('are your points in a plane?')

x, y, z = np.indices((xdim, ydim, zdim))

# use the plane equation dot(n,x) = c to find where the 
plane = abs(n[0] * x + n[1] * y + n[2] * z - c) <= 0.5

# color the plane blue ...
facecolors = np.empty(plane.shape, dtype=object)
facecolors[plane] = 'blue'

# ... and the given points red
facecolors[x0, y0, z0] = 'red'
facecolors[x1, y1, z1] = 'red'
facecolors[x2, y2, z2] = 'red'
# make sure they will be filled, just in case
plane[x0, y0, z0] = True
plane[x1, y1, z1] = True
plane[x2, y2, z2] = True

fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.voxels(plane, facecolors=facecolors);
plt.show()

Screenshot, for a case with 3 vertices visible:

3D hyper-plane in a numpy array

1 Answers1