Given a plane through the origin with normal \(\mathbf{n}\) and a reflect direction \(\mathbf{r}\), we are looking for a transformation matrix \(\mathbf{T}\) that is able to calculate the reflected ray from the plane such that \(\mathbf{r}' = \mathbf{T}\cdot\mathbf{r}\)
The incoming direction \(\mathbf{r}\) can be expressed as two components, a vector \(\mathbf{r}^|\) parallel to the normal vector \(\mathbf{n}\) and one \(\mathbf{r}^\perp\), perpendicular to it:
\[\mathbf{r} = \mathbf{r}^| + \mathbf{r}^\perp\]
The parallel vector is simply the vector projection of \(\mathbf{r}\) onto \(\mathbf{n}\):
\[\mathbf{r}^| = \text{proj}_\mathbf{n}(\mathbf{r})=\frac{\mathbf{r}\cdot\mathbf{n}}{\mathbf{n}\cdot\mathbf{n}}\mathbf{n}\]
And the perpendicular is simply
\[\mathbf{r}^\perp = \mathbf{r}-\mathbf{r}^|\]
When we now reflect the vector \(\mathbf{r}\), the component parallel to the normal remains unchanged, while the component perpendicular to it is negated:
\[\mathbf{r}^{'\perp} = -\mathbf{r}^\perp\]
Now combining this information leads to
\[\mathbf{r}' = \mathbf{r}^| + \mathbf{r}^{'\perp} = \mathbf{r}^| - \mathbf{r}^\perp\]
Plugging in and simplifying leads to
\[\begin{array}{rl} \mathbf{r}' &= \mathbf{r}^| - \mathbf{r}^\perp\\ &= \frac{\mathbf{r}\cdot\mathbf{n}}{\mathbf{n}\cdot\mathbf{n}}\mathbf{n} - \left(\mathbf{r} - \frac{\mathbf{r}\cdot\mathbf{n}}{\mathbf{n}\cdot\mathbf{n}}\mathbf{n}\right)\\ &= \mathbf{r} - 2\frac{\mathbf{r}\cdot\mathbf{n}}{\mathbf{n}\cdot\mathbf{n}}\mathbf{n}\\ \end{array}\]
Since we are looking for a transformation matrix \(\mathbf{T}\) to express the reflection, we simplify the equation by normalizing the normal \(\hat{\mathbf{n}}\) and factor out \(\mathbf{r}\) to get the desired matrix \(\mathbf{T}\):
\[\mathbf{r}' = \underbrace{(\mathbf{I} - 2\hat{\mathbf{n}}\cdot\hat{\mathbf{n}}^T)}_{=\mathbf{T}}\mathbf{r}\]
Reflection on Normal
When looking at the equation \(\mathbf{r}'=\mathbf{r} - 2\frac{\mathbf{r}\cdot\mathbf{n}}{\mathbf{n}\cdot\mathbf{n}}\mathbf{n}\) again, it is interesting that this basically says that
\[\begin{array}{rl} \mathbf{r}' &=\mathbf{r} - 2\frac{\mathbf{r}\cdot\mathbf{n}}{\mathbf{n}\cdot\mathbf{n}}\mathbf{n}\\ & = \mathbf{r} - 2\cdot\text{proj}_\mathbf{n}(\mathbf{r})\\ & = -(2\cdot\text{proj}_\mathbf{n}(\mathbf{r}) - \mathbf{r})\\ & = -\text{refl}_\mathbf{n}(\mathbf{r}) \end{array}\]
or in other words that our reflected vector \(\mathbf{r}'\) for the plane is the negated vector of what is typically called the reflection vector \(\text{refl}_\mathbf{n}(\mathbf{r})\) on vector \(\mathbf{n}\).