Get list of vertices from List of planes

517 Views Asked by At

Hello I have an example 3d object defined by these seven planes. By my reckoning, all solids produced in this way will be convex (if it matters).
( (-512 0 512) (-512 256 512) (0 256 512) )
( (-512 512 0) (-512 0 0) (0 0 0) )
( (-512 0 0) (-512 512 0) (-512 512 256) )
( (0 512 0) (0 0 0) (0 0 512) )
( (-512 512 0) (0 512 0) (0 512 256) )
( (0 0 0) (-512 0 0) (-512 0 512) )
( (0 512 256) (0 256 512) (-512 256 512) )

Is there a way to get a list of the vertices formed by the intersection of these planes, removing the two vertices created where the top plane meets the front right plane in the image below.

2D Simplification

This pentagon is dedfined by the following 5 lines y = 0, y = 2, x = 0, x = 2 and y = -x + 3

Is there a way to return points A-E, but Not F from the 5 lines, and is that method scalable to 3 dimensions?

Many thanks

1

There are 1 best solutions below

4
On BEST ANSWER

The planes you refer to are called half-spaces. A polytope defined by the intersection of half-spaces is indeed always convex.

First, convert your point triplets (defining a plane) to the form with a normal vector $\vec{n}_i$ pointing outwards from the polytope, and signed distance $d_i$ from origin. In this form, point $\vec{p}$ is inside the half-space, if $$\vec{p} \cdot \vec{n}_i \le d_i$$ and point $\vec{p}$ is inside the polytope, if and only if the above holds for all half-spaces $i$.

(Note that in this notation, $d_i$ is negated compared to what e.g. the Wikipedia page uses for example for $a x + b y + c z + d \le 0$.)


Calculating $\vec{n}_i$ and $d_i$, when three non-collinear points $\vec{p}_1 = (x_1, y_1, z_1)$, $\vec{p}_2 = (x_2, y_2, z_2)$, and $\vec{p}_3 = (x_3, y_3, z_3)$ are known, is quite straightforward.

The Wikipedia Plane (geometry) describes an easy way to obtain $\vec{n}_i$: $$\vec{n}_i = \left(\vec{p}_2 - \vec{p}_1\right) \times \left(\vec{p}_3 - \vec{p}_1\right)$$ When $\vec{n}_i$ is known, finding out $d_i$ is trivial: $$d_i = \vec{p}_1 \cdot \vec{n}_i$$ Above, you can use any point on the plane instead of $\vec{p}_1$, because $d_i$ is just the signed distance from the plane to origin (at minimum magnitude, i.e. closest point).

Note that you also need to verify that each $\vec{n}_i$ points away from the object. If the three points specifying each plane are defined in counterclockwise order if viewed from outside, then the vectors should all point outwards. Otherwise, you need to check by hand. To switch the direction of the normal vector, negate $d_i$ and each component of $\vec{n}_i$.


To find all vertices, you need to find the intersection of each triplet of non-parallel half-spaces. If the intersection point is inside the polytope, it is a vertex. (If you also need the edge list, you should include the indexes $i$ of the half-spaces that defined the vertex -- and you do need to combine duplicate vertices, but keep all the unique indices -- as the set of vertices that have index $i$ define the face $i$; you only need to order them properly to determine the exact edges.)

If we have $\vec{n}_i = (x_i, y_i, z_i)$, $\vec{n}_j = (x_j, y_j, z_j)$, and $\vec{n}_k = (x_k, y_k, z_k)$, then we can calculate their intersection $(x, y, z)$ quite easily. First, we need to calculate a temporary: $$\begin{align}t &= x_i y_j z_k - x_i z_j y_k - y_i x_j z_k + y_i z_j x_k + z_i x_j y_k - z_i y_j x_k \\ &= x_i ( y_j z_k - z_j y_k ) + y_i ( z_j x_k - x_j z_k ) + z_i ( x_j y_k - y_j x_k )\end{align}$$ If the planes are parallel, then $t = 0$. Otherwise, $$\vec{p} = \begin{cases} x = \frac{d_i (z_j y_k - y_j z_k) + d_j (y_i z_k - z_i y_k) + d_k (z_i y_j - y_i z_j)}{t} \\ y = \frac{d_i (x_j z_k - z_j x_k) + d_j (z_i x_k - x_i z_k) + d_k (x_i z_j - z_i x_j)}{t} \\ z = \frac{d_i (y_j x_k - x_j y_k) + d_j (x_i y_k - y_i x_k) + d_k (y_i x_j - x_i y_j)}{t} \end{cases}$$ If and only if $\vec{p} \cdot \vec{n}_m \le d_m$ for all $m$ is $\vec{p}$ a vertex. This final check excludes points such as $F$ in OP's diagram.

In pseudocode, for $N$ half-spaces:

eps = 0.000001

For i = 1 to (N-2):
    For j = (i+1) to (N-1):
        For k = (j+1) to N:
            t = x[i]*(y[j]*z[k] - z[j]*y[k]) + y[i]*(z[j]*x[k]-x[j]*z[k]) + z[i]*(x[j]*y[k]-y[j]*x[k])
            If (abs(t) <= eps):
                continue
            End if

            x = (d[i]*(z[j]*y[k]-y[j]*z[k]) + d[j]*(y[i]*z[k]-z[i]*y[k]) + d[k]*(z[i]*y[j]-y[i]*z[j])) / t
            y = (d[i]*(x[j]*z[k]-z[j]*x[k]) + d[j]*(z[i]*x[k]-x[i]*z[k]) + d[k]*(x[i]*z[j]-z[i]*x[j])) / t
            z = (d[i]*(y[j]*x[k]-x[j]*y[k]) + d[j]*(x[i]*y[k]-y[i]*x[k]) + d[k]*(y[i]*x[j]-x[i]*y[j])) / t

            For m = 1 .. N:
                If (m != i) and (m != j) and (m != k) and
                   (x*x[m] + y*y[m] + z*z[m] > d[m] + eps):
                    m = -1
                    Break
                End if
            End for
            If (m >= 1):
                # Point (x,y,z) is a vertex.
            End if
        End for
    End for 
End for

The above is essentially $O(N^4)$ and thus quite slow, but for smallish number of faces it should do just fine. (Actually, the inner math is fast enough for it to be surprisingly usable for quite complex convex objects.)

The use of eps is a reminder that floating-point values are not exact, and there is always some rounding error; instead of comparing to exact values, one should compare the magnitude to some epsilon.