Given a plane through the origin with normal \(\mathbf{n}\) and an incoming ray direction \(\mathbf{r}\), we are looking for a transformation matrix \(\mathbf{T}\) that gives the reflected direction directly as \(\mathbf{r}' = \mathbf{T}\cdot\mathbf{r}\).
Decomposing the Ray into Parallel and Perpendicular Parts
The incoming direction \(\mathbf{r}\) can be split into two components, a vector \(\mathbf{r}_\parallel\) parallel to the normal vector \(\mathbf{n}\) and one \(\mathbf{r}_\perp\), perpendicular to it, i.e. lying in the plane itself:
\[\mathbf{r} = \mathbf{r}_\parallel + \mathbf{r}_\perp\]
The parallel part is simply the vector projection of \(\mathbf{r}\) onto \(\mathbf{n}\):
\[\mathbf{r}_\parallel = \text{proj}_\mathbf{n}(\mathbf{r}) = \frac{\mathbf{r}\cdot\mathbf{n}}{\mathbf{n}\cdot\mathbf{n}}\mathbf{n}\]
and the perpendicular part is what's left over:
\[\mathbf{r}_\perp = \mathbf{r} - \mathbf{r}_\parallel\]
Reflecting Across the Plane
Reflecting \(\mathbf{r}\) across the plane leaves everything that already lies in the plane untouched and only flips the part sticking out of it. That is exactly \(\mathbf{r}_\perp\) and \(\mathbf{r}_\parallel\) respectively: the in-plane component is unchanged, while the component along the normal is negated,
\[\mathbf{r}'_\perp = \mathbf{r}_\perp,\qquad \mathbf{r}'_\parallel = -\mathbf{r}_\parallel.\]
Adding both reflected components back together gives the reflected ray
\[\mathbf{r}' = \mathbf{r}'_\parallel + \mathbf{r}'_\perp = \mathbf{r}_\perp - \mathbf{r}_\parallel.\]
Plugging in \(\mathbf{r}_\perp = \mathbf{r} - \mathbf{r}_\parallel\) and simplifying leads to
\[\begin{array}{rl} \mathbf{r}' &= \mathbf{r}_\perp - \mathbf{r}_\parallel\\ &= \left(\mathbf{r} - \frac{\mathbf{r}\cdot\mathbf{n}}{\mathbf{n}\cdot\mathbf{n}}\mathbf{n}\right) - \frac{\mathbf{r}\cdot\mathbf{n}}{\mathbf{n}\cdot\mathbf{n}}\mathbf{n}\\ &= \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 normalize the normal to \(\hat{\mathbf{n}}\) and factor \(\mathbf{r}\) out of the expression 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}\]
Relation to Reflecting on a Vector
Looking at \(\mathbf{r}' = \mathbf{r} - 2\frac{\mathbf{r}\cdot\mathbf{n}}{\mathbf{n}\cdot\mathbf{n}}\mathbf{n}\) again, it is interesting that this is basically saying
\[\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, 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 the vector \(\mathbf{n}\). Reflecting on a plane and reflecting on its normal produce the same in-plane behavior but opposite outcomes along the normal — which is exactly why one is the negation of the other.
Implementation
Turning \(\mathbf{r}' = \mathbf{r} - 2\frac{\mathbf{r}\cdot\mathbf{n}}{\mathbf{n}\cdot\mathbf{n}}\mathbf{n}\) into code works component-wise for any dimension, without ever assembling the matrix \(\mathbf{T}\) explicitly:
function reflectOnPlane(r, n) {
const dot = (a, b) => a.reduce((sum, ai, i) => sum + ai * b[i], 0);
const k = 2 * dot(r, n) / dot(n, n);
return r.map((ri, i) => ri - k * n[i]);
} Reflecting \(\mathbf{r}=(-3,3,0)\) on a plane with normal \(\mathbf{n}=(0,1,0)\), for instance, gives \(\mathbf{r}'=(-3,-3,0)\): the component along the normal flips sign, while the component in the plane is left alone.