The Orthogonal Basis
The first orthorgonal vector to $\mathbf{V}_1$, called $\mathbf{V}_1$ is done by
$$
\mathbf{V}_2=(-V_y, V_x, 0).
$$
We can easily verify that $\mathbf{V}_1$ and $\mathbf{V}_2$ are orthogonal by taking the inner product
$$
\langle\mathbf{V}_1,\mathbf{V}_2\rangle = V_x\cdot(-V_y)+V_y\cdot V_x + V_z\cdot 0 = -V_xV_y+V_xV_y = 0.
$$
The last orthogonal vector $\mathbf{V}_3$ is created by taking the cross-product of $\mathbf{V}_1$ and $\mathbf{V}_2$ yielding
$$
\mathbf{V}_3=\mathbf{V}_1\times\mathbf{V}_2=(V_y\cdot 0 - V_z\cdot V_x,V_z\cdot (-V_y)-V_x\cdot 0,V_x\cdot V_x - V_y\cdot (-V_y)),
$$
and can be reduced to
$$
\mathbf{V}_3=(-V_x\cdot V_z, -V_y\cdot V_z, V^2_x+V^2_y).
$$
Thus $\mathbf{V}_1$, $\mathbf{V}_2$ and $\mathbf{V}_3$ now constitue an orthogonal basis. However, normalizing by brute-force is cumbersome (and error prone) so we will take a slightly more clever route.
The Orthonormal Basis
If we instead assume that the vector $\mathbf{V}_1$ is normalized initially, we have that$$
\hat{\mathbf{V}}_1 = (V_x, V_y, V_z).
$$
I've re-used the same symbols as before to reduce an unnecessary cluttering of the notation, but do remember that elements of $\hat{\mathbf{V}}_1$ are normalized. To help us later on, we can write the length of $\hat{\mathbf{V}}_1$ as
$$
1=\sqrt{V^2_x + V^2_y + V^2_z}=> 1=V^2_x + V^2_y + V^2_z
$$
which we will use momentarily. $\hat{\mathbf{V}}_2$ is constructed as above but we have to normalize right away
$$
\hat{\mathbf{V}}_2 = (\frac{-V_y}{\sqrt{V^2_x + V^2_y}}, \frac{V_x}{\sqrt{V^2_x + V^2_y}}, 0)=(\frac{-V_y}{\sqrt{1-V^2_z}}, \frac{V_x}{\sqrt{1-V^2_z}}, 0)
$$
In the final equality sign, we used the length of $\hat{\mathbf{V}}_1$.
Now, for the cross-product. Remember that the no normalization for $\hat{\mathbf{V}}_3$ is needed because $\hat{\mathbf{V}}_1$ and $\hat{\mathbf{V}}_2$ are already normalized
$$
\hat{\mathbf{V}}_3=\hat{\mathbf{V}}_1\times\hat{\mathbf{V}}_2=(\frac{-V_z\cdot V_x}{\sqrt{1-V^2_z}},\frac{-V_z\cdot V_y}{\sqrt{1-V^2_z}},\frac{V^2_x+V^2_y}{\sqrt{1-V^2_z}}).
$$
The last element of the $\hat{\mathbf{V}}_3$ vector is something we can take care of and make it a tad prettier. By using the length of vector $\hat{\mathbf{V}}_1$ again we obtain
$$
\hat{\mathbf{V}}_3=(\frac{-V_z\cdot V_x}{\sqrt{1-V^2_z}},\frac{-V_z\cdot V_y}{\sqrt{1-V^2_z}},\frac{1-V^2_z}{\sqrt{1-V^2_z}}),
$$
and if we remember that $x/\sqrt(x)=\sqrt(x)$ we finally obtain $\hat{\mathbf{V}}_3$ as
$$
\hat{\mathbf{V}}_3=(\frac{-V_z\cdot V_x}{\sqrt{1-V^2_z}},\frac{-V_z\cdot V_y}{\sqrt{1-V^2_z}},\sqrt{1-V^2_z}).
$$
So to sum up, given a normalized vector $\hat{\mathbf{V}}_1$, we can construct an orthonormal basis using the following equations
$$
\hat{\mathbf{V}}_1 = (V_x, V_y, V_z),
$$
$$
\hat{\mathbf{V}}_2 =(\frac{-V_y}{\sqrt{1-V^2_z}}, \frac{V_x}{\sqrt{1-V^2_z}}, 0),
$$
$$
\hat{\mathbf{V}}_3=(\frac{-V_z\cdot V_x}{\sqrt{1-V^2_z}},\frac{-V_z\cdot V_y}{\sqrt{1-V^2_z}},\sqrt{1-V^2_z}).
$$
The clever thing about these last expressions is that the normalization factor $\sqrt{1-V^2_z}$ is present many times and can be precomputed for each new $\hat{\mathbf{V}}_1$.
Acknowledgements: Thanks to Lars for being awesome with the math.
