Find extra arbitrary two points for a plane, given the normal and a point that lies on the plane

Question

For a plane, I have the normal $n$, and also a point $P$ that lies on the plane.

Now, how am I going to find extra arbitrary two points ($P_1$ and $P_2$) for the plane so that these three points $P$, $P_1$ and $P_2$ completely define the plane?

The solution here suggests that one assumes a certain $x$ and $y$ to substitute into the plane equation and find the remaining $z$. But this method is only suitable for hand calculation; it breaks down for plane $z=0$. As such, it is not suitable for computer implementation.

I would need an algorithm that is robust and can handle all the cases, any idea how to construct it?

There is a similar question here, and the answer suggests me to use Gram-Schmidt, but I don't see how it can be applied in my case here.

Answer 1

The fundamental problem here is to find two vectors $u, v$ that are orthogonal to the vector $n = (n_x, n_y, n_z)^T$. Morevover, the procedure should be numerically stable. I assume that $|n| = 1$.

Approach 1. Find the component of $n$ that has the smallest absolute value, say $n_z = \epsilon_z$. The vector $n$ is hence "almost" in the x/y-plane (z = 0). Take the vector $w = (0, 0, 1)^T$, with all components zero except that one on the place of the smallest component. Compute $u = w \times n$. $u$ is orthogonal to $n$ and is computed numerically stable because it is "almost" orthogonal to $n$. Finalize the procedure by computing $v = n \times u$.

Approach 2. Often Householder orthogonalization is recommended for ill-conditioned orthogonalization. It could work better than the Approach 1.

1. Build that "mirror" vector $w = (n_x + 1, n_y, n_z)$ or $w = (n_x - 1, n_y, n_z)$ whose first component has bigger absolute value.
2. Compute the 3x3 Householder matrix $H = I - 2 w w^T / (w^T w)$
3. The first row of H will be a unit vector parallel to $n$, and the other two rows will be unit vectors orthogonal to $n$ and to each other, i.e. these are the vectors $u^T, v^T$.

After $u, v$ are computed, we set $P_1 = P + u, P_2 = P + v$.

Answer 2

Case 1: If the plane does not pass the origin, then there is no basis of this plane because the plane does not form a subspace.

Case 2: If the plane passes the origin, then only one normal vector $n$ is already sufficient to determine this plane, and there is no need to give the point $P$.

In the case 2, let $x$ be an point on the plane, then $$n^Tx=0$$ Using Singular Value Decomposition (SVD), you can find an orthogonal basis of the plane. The procedure is as below. Suppose one SVD of $n\in\mathbb{R}^d$ is $$n=U\Sigma V^T$$ where $U\in\mathbb{R}^{d\times d}$ and $V\in\mathbb{R}$. Denote $U=(u_1,u_2,\cdots,u_d)$, then $$u_i^Tn=n^Tu_i=0$$ for all $i=2,\cdots,d$. So one orthogonal basis of the plane is $\{u_2,\cdots,u_d\}$.

Answer 3

Let $n=(a,b,c)$. If $a=b=0$, take $(1,0,0)$ and $(0,1,0)$ as the basis. Otherwise, take $u=(b,-a,0)$ and $v=n \times u$.

Answer 4

Here is a naive "hack" that will likely work, and will avoid complex computations. My apologies to the purists.

First, that the plane may not pass through the origin is irrelevant, as others have mentioned. Just imagine moving it so it does pass through the origin (e.g., by subtracting $p$ from every point on the plane), and after computing the $p_1$ and $p_2$ you want for the shifted plane, add $p$ to both points to shift the plane back.

Second, you only need find $p_1$ robustly, because once you have $p_1$ not too close to the origin (the former $p$), you can just rotate $p_1$ about $n$ by $90^\circ$.

So it all comes down to finding a point $p_1$ on a plane through the origin that is not itself too close to the origin, given the normal $n$. Let $n=(a,b,c)$ as per lhf. Any point $(x,y,z)$ on the plane satisfies $a x + b y + c z = 0$.

Here's the hack. Choose a random point on the plane, say by selecting two coordinates, generating random numbers, and solving for the third coordinate. Of course you have to be a bit intelligent about the cases where one or two coordinates are zero, but that is not too difficult. If you end up too close to the origin for your robustness requirements, try again.

Eventually you find $p_1$. Rotate about $n$ to get $p_2$. Finally shift back to $\{ p, p_1 + p, p_2 + p \}$.