Find extra arbitrary two points for a plane, given the normal and a point that lies on the plane 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. 
 A: 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.


*

*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.

*Compute the 3x3 Householder matrix $H = I - 2 w w^T / (w^T w)$

*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$.
A: 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\}$.
A: 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$.
A: 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 \}$.
A: If you any component of your normal is non-zero (and we know at least ONE component won't be), any value of the other two components will fall SOMEWHERE on the plane.  Say you have a normal with a positive or negative y component, you know any pair of x, z values will fall somewhere on its plane.  Plug in any two values for it and solve for the last using the plane equation.
vec3 get_point_on_plane(vec3 const& point, vec3 const& norm)
{
//get random x, y, z points
float x = point.x + 1;
float y = point.y + 1;
float z = point.z + 1;

/*
    Plane equation:
    Ax + By + Cz = D  ==  
    norm.x(x - point.x) + norm.y(y - point.y) + norm.z(z - point.z) = 0
*/
if (! nearly_equal(norm.x, 0)) //any y, z pair will fall SOMEWHERE on the plane
{
    //solve for x with rand_y and rand_z
    x = point.x - ( norm.y * (y - point.y) + norm.z * (z - point.z) ) / norm.x;

}
else if (! nearly_equal(norm.y, 0)) //any x, z pair will fall SOMEWHERE on the plane
{
    //solve for y with rand_x and rand_z
    y = point.y - ( norm.x * (x - point.x) + norm.z * (z - point.z) ) / norm.y;
}
else //any x, y pair will fall SOMEWHERE on the plane
{
    //solve for z with rand_x and rand_y
    z = point.z - ( norm.x * (x - point.x) + norm.y * (y - point.y) ) / norm.z;
}
return vec3(x, y, z);

}
