Check Whether a 2D Plane Intersects a Hypercube

Question

I'd like to quickly check whether a 2D plane intersects a $N$-dimensional hypercube. In my case, the hypercube is $[0,1]^N$, and the plane is described by an offset point $\mathbf r$ and two vectors $\mathbf u$ and $\mathbf v$.

$$\mathbf r = \mathbf r_0 + s \mathbf u + t \mathbf v$$

I'm not interested in constructing the intersection or identifying any of its properties, only determining whether it exists.

I didn't find an answer in these related questions:

• This question concerns the intersection of a random subspace centered on the origin with a hypercube. My application involves a specified plane that may not intersect the origin.
• These questions 1, 2 concern properties of the intersection, but do not touch upon fast methods for testing whether the intersection is nonempty.
• This question concerns checking if a line intersects a hypercube, but I wasn't quite sure how to generalize it for checking if a plane intersects a hypercube.

Edit:

I've removed a more extended example discussing a possible answer which has since been deleted.

The current answer rightly suggests to use Linear Programming. However, I'd hoped that there was a more direct solution based on the geometry of the problem, something vaguely akin to this approach for testing whether a point lies inside a triangle.

Answer 1

This can be formulated as a linear programming problem. Break the problem up by coordinates.

\begin{align*} 0 &\leq (\mathbf{r}_0 + s \mathbf{u} + t \mathbf{v})^{1} \leq 1 \\ 0 &\leq (\mathbf{r}_0 + s \mathbf{u} + t \mathbf{v})^{2} \leq 1 \\ & \phantom{(r_0+s u+)}\vdots \\ 0 &\leq (\mathbf{r}_0 + s \mathbf{u} + t \mathbf{v})^{N} \leq 1 \\ \end{align*}

Let $\mathbf{x} = (s,t)^{\text{T}} \in \Bbb{R}^2$. You only want to know if there is a feasible point satisfying these inequalities; you aren't trying to maximize or minimize anything. So pick your favorite or a random objective function "to maximize". Then convert each inequality, for $i=1\dots N$, above into the pair of feasibility constraints \begin{align*} 0 &\leq (\mathbf{r}_0 + s \mathbf{u} + t \mathbf{v})^{i} \\ -(\mathbf{r}_0 + s \mathbf{u} + t \mathbf{v})^{i} &\leq 0 \\ -\mathbf{r}_0^i - s \mathbf{u}^i - t \mathbf{v}^{i} &\leq 0 \\ - s \mathbf{u}^i - t \mathbf{v}^{i} &\leq \mathbf{r}_0^i \end{align*} and \begin{align*} (\mathbf{r}_0 + s \mathbf{u} + t \mathbf{v})^{i} &\leq 1 \\ \mathbf{r}_0^i + s \mathbf{u}^i + t \mathbf{v}^{i} &\leq 1 \\ s \mathbf{u}^i + t \mathbf{v}^{i} &\leq 1 -\mathbf{r}_0^i \text{.} \end{align*} So the coefficient matrix, $A$, has $2N$ rows in the pairs \begin{pmatrix} &\vdots \\ -\mathbf{u}^i & &-\mathbf{v}^i \\ \mathbf{u}^i & & \mathbf{v}^i \\ &\vdots \end{pmatrix} and the corresponding rows of $\mathbf{b}$ are \begin{pmatrix} \vdots \\ \mathbf{r}_0^i \\ 1 - \mathbf{r}_0^i \\ \vdots \end{pmatrix} This makes the constraints equation $A \mathbf{x} \leq \mathbf{b}$.

Linear programming makes the additional constraints $\mathbf{x} \geq 0$, which corresponds to finding an intersection in one quadrant of the given plane. Since we want to check all four quadrants (but are free to stop once we find any point at all), we run an LP solver up to four times. The first, as described above. The second, with $\mathbf{u}$ replaced with $-\mathbf{u}$; this is to implement the replacement $s \mapsto -s$ to search for solutions with $s \leq 0$ coordinate. Then with $\mathbf{v}$ replaced with $-\mathbf{v}$, to search for solutions with $t \leq 0$. Then a fourth time with both $\mathbf{u}$ and $\mathbf{v}$ replaced with their negatives to search the remaining quadrant.

If any run finds a feasible point, there is an intersection. As soon as a run finds a feasible point, you need no further runs. If all four runs find no feasible points, then there is no intersection.

It turns out that finding a feasible point is about as hard (equivalent computational complexity) as solving the LP instance. Some LP solvers will let you stop between various phases of their computation. If you use a solver that will let you stop as soon as it finds a feasible point (any point in the intersection of the cube and (quadrant of the) plane), you can save some run time.

(I spent a little while trying to leverage the $x \leq 0$ constraints to be half of the constraints for the cube, so we would not need up to four LP runs. This would have $\mathbf{x} \in \Bbb{R}^N$. The obstruction is that there does not appear to me to be a linear inequalities way to determine whether a particular choice of $\mathbf{x}$ is on the plane.)

Answer 2

The answer based on Linear Programming (LP) is is clear and easy to implement, and can be used generally for a much broader class of intersection problems. However, we can gain a bit of efficiency in some cases by exploiting the geometry of this problem directly, as described below.

Method 1: direct test for intersection

This is fast for low dimensions e.g. N=4 or 5, but scales poorly for large N.

Consider projecting a unit $N$-cube centered at $\mathbf p \in \mathbb R^N$ onto a plane centered at the origin and defined by the vectors $\mathbf u, \mathbf v \in \mathbb R^N$. (The more general problem of checking for the intersection of a $N$-cube and a plane can be reduced to this with the right choice of coordinates.)

The plane defines a 2D $(s,t)$ subspace of $N$-D space, with basis $A = (\mathbf u, \mathbf v)$

$$ \mathbf x = A \begin{bmatrix} s\\t \end{bmatrix} $$

We can project the problem onto the null space of the plane, $A^\perp$. This sends every point on the plane to zero. The problem then reduces to testing whether the $N{-}2$ projection of the hypercube contains the origin.

Projecting the $N$-cube down to $N{-}2$ dimensions creates a polytope that can be expressed as the union of all $N{-}2$ faces of the original $N$-cube. If any of these faces contains the origin, then the original hypercube intersects the plane.

E.g. projecting a 4-cube down to 2D yields a collection of 2D rhombi, one for each 2D faces of the 3D facets of the 4-cube. If any of these rhombi contain the origin, then we know the 2D plane intersects the 4-cube in our original problem.

One way to check if a $N{-}2$ rhombus contains the origin is to change to a basis where the rhombus is a cube $[0,1]^{N-2}$, and test if this cube contains the origin. Do this by selecting one vertex of the rhombus, and all the points it immediately connects to, as your basis set.

(When checking all sub-facets, you can stop as soon as you've found a single sub-facet that contains the target point)

This seems inefficient, but is faster than using e.g. Python's built-in linear programming solvers for $N{=}4$. It scales badly to higher dimensions, however. There might be a more elegant solution if one could exploit symmetries and eliminate redundant computations.

Method 2: $\mathcal O[ N \log(N) ]$ via Preparata and Muller's algorithm (or something like it)

First, transform the problem into a set of $2N$ linear inequality constraints. Each constraint defines a half-plane. The plane intersects the hypercube if the intersection of these half-planes is nonempty.

So far, this is identical to the linear programming approach. But, as it turns out, there are specific algorithms for testing whether the intersection of $n=2N$ half-planes is nonempty. For example, here is an approach by Preparata and Muller for solving the intersection of $n$ half-spaces with $n\log(n)$ time complexity. There are other variants (e.g. Wu, Ji, and Chen), but they all have the same complexity.

These lecture notes by Dave Mount are especially useful for understanding the geometry underlying these algorithms. These notes show how to construct the (convex) intersection set as the intersection of a (convex) upper and lower envelope.

The basic pseudocode is:

• First, identify any vertical bounding lines. These demarcate spans of the $s$ axis of the plane. Their intersection $s\in[s_0,s_1]$ defines the bounds for a search procedure (below).
• Split the remaining lines into those that bound the half-plane from below, and those that bound the half plane from above, where "below" and "above" are defined in terms of the $t$ coordinate of the plane.
• These two sets of lines define and upper and lower feasibility regions, which are convex. The boundary of these regions can be interpreted as curves $t_l(s)$ and $t_u(s)$
• We can test if the the intersection between the upper/lower feasibility regions is nonempty by finding the minimum of $t(s) = t_l(s)-t_u(s)$ on $s\in[s_0,s_1]$
• If any $\exists s\in[s_0,s_1]\text{ s.t. } t(s)<0 $, then the plane intersects the hypercube.
• This can be checked via binary search, looking for the point where $t(s)$ changes sign, and stopping early if any point satisfying all the constraints is found.