The answer in the screenshot is wrong on two counts. First, you're correct that in the middle matrix the −1 should be on the third row, not the second. The other error is that the $u, v$ basis vectors are normalized, but not orthogonalized, so $M$ is not in general a rotation matrix; therefore its inverse is not simply its transpose. (The reflection doesn't actually need an orthonormal basis, so you could fix the answer either way: by using the full inverse for $M$, or by fully orthonormalizing $u$ and $v$.)
The rule of thumb for setting up a local coordinate system (I assume orthonormal) for a given object is to set some of the axes to be lined up with the object, then choose arbitrary perpendicular vectors for any remaining axes.
- For a point, there's nothing to line up with, so you can use any arbitrary axes.
- For a line, you line up the $u$ axis with the line, then there's nothing to constrain the $v$ axis, so you can choose any arbitrary vector perpendicular to $u$. That's what the $v=u\times(0,0,1)^T$ thing was getting at: it's a convenient formula to get some vector perpendicular to $u$ if you don't otherwise care what direction it is. Except that this won't work if $u$ happens to be parallel to the z-axis (the cross product returns zero then), so you actually need to detect that and switch to a different formula, such as $v=u\times(1,0,0)^T$, in that case.
- For a plane, you line up the $u$ and $v$ axes with the plane. So, you can pick any arbitrary $u$ that lies in the plane, then find a $v$ that lies in the plane and is perpendicular to $u$.
Addendum on orthonormalization
Suppose you have a couple of 3D vectors $u, v$ that are non-parallel (linearly independent) but otherwise arbitrary, and you want to generate a 3D orthonormal basis $u', v', w'$ from them. There are two equivalent ways to do it.
Use the Gram–Schmidt process. This works as follows: first, normalize $u$. Then, for each additional vector, subtract off its projection onto the previous vectors, and normalize what's left over. For the 3D case this looks like:
$$ \begin{aligned}
u' &= \text{normalize}(u) = \frac{u}{|u|} \\
v' &= \text{normalize} \bigl( v - (u' \cdot v)u' \bigr) \\
w' &= u' \times v'
\end{aligned} $$
Note that the notation $(u' \cdot v)u'$ means take the dot product of $u'$ with $v$ (giving a scalar), then multiply that scalar by $u'$. This gives the vector component of $v$ that is parallel to $u'$. When this is subtracted from $v$, the remaining component is perpendicular to $u'$.
Finally, $w'$ is set to the cross product of $u'$ and $v'$ to ensure the basis is right-handed. Since $u', v'$ are already orthonormal, we don't need to normalize $w'$—it comes out already normalized. (This isn't strictly part of the Gram–Schmidt process.)
By the way, the Gram–Schmidt process works for making an orthonormal basis in any number of dimensions, not just two or three.
The other way, which only works in 3D, is to use cross products.
$$ \begin{aligned}
u' &= \text{normalize}(u) \\
w' &= \text{normalize}(u' \times v) \\
v' &= w' \times u'
\end{aligned} $$
Here, $w'$ has to be normalized because $u'$ and $v$ are not orthonormal. But then $u'$ and $w'$ are orthonormal, so we can calculate $v'$ from them without a normalize. This gives the same results as the previous process, so it's up to you which you prefer to use.
