3D reflection matrix

3d_reflection_matrix.txt

Affine transformation matrices to reflect a point in (1) the XY-plane through the origin; (2) an arbitrary 3D plane. 
A right-handed coordinate system is assumed.

1. XY plane at z=0

Consider a reflection in an XY-plane located at the origin. To accomplish this we require a matrix M such that

  |x|   | x|
M |y| = | y|
  |z|   |-z|
  |1|   | 1|

Such an affine matrix is

    | 1   0   0   0 |
M = | 0   1   0   0 |
    | 0   0  -1   0 |
    | 0   0   0   1 |



2. Arbitrary 3D plane

Now the plane and point it reflects will first be subject to a sequence of reversible transformations such that after they
are applied the plane coincides with the origin and is orientated so that its normal vector points in the positive Z-axis 
direction. Since the plane now coincides with the XY plane at z=0 the actual reflection transformation which reflects the 
point in the plane is applied by using the XY-plane reflection matrix. Finally, the transformations which were applied to 
change the plane's position and orientation are reversed.

See https://gist.github.com/fospathi/1bfbc5d5a1f1df588aa31732a8cd451e for more details.

Let the plane be defined by a point on the plane Q=(q_x, q_y, q_z) and a normal direction D=(d_x, d_y, d_z).

Let R be a transformation matrix which reflects in the plane, then we have

    |1 0 0 q_x| | cz sz 0 0| |cy 0 -sy 0| |1 0  0 0| | cy 0 sy 0| |cz -sz 0 0| |1 0 0 -q_x|
R = |0 1 0 q_y|*|-sz cz 0 0|*|0  1  0  0|*|0 1  0 0|*| 0  1 0  0|*|sz  cz 0 0|*|0 1 0 -q_y|
    |0 0 1 q_z| | 0  0  1 0| |sy 0  cy 0| |0 0 -1 0| |-sy 0 cy 0| |0   0  1 0| |0 0 1 -q_z|
    |0 0 0 1  | | 0  0  0 1| |0  0  0  1| |0 0  0 1| | 0  0 0  1| |0   0  0 1| |0 0 0  1  |

Let M be the product of the leftmost six of these seven matrices, then

        |1 0 0 -q_x|
R = M * |0 1 0 -q_y|
        |0 0 1 -q_z|
        |0 0 0  1  |

where 

θ_y = -atan2{√[(d_x)²+(d_y)²], d_z}
θ_z = -atan2{d_y, d_x}

cy = cos{θ_y}
cz = cos{θ_z} 
sy = sin{θ_y}
sz = sin{θ_z}

M₀₀ = cy*cy*cz*cz + sz*sz - cz*cz*sy*sy
M₀₁ = cz*sz*(1-cy*cy) + cz*sy*sy*sz
M₀₂ = 2 * cy * cz * sy
M₀₃ =  q_x
M₁₀ = cz*sz*(1-cy*cy) + cz*sy*sy*sz
M₁₁ = cy*cy*sz*sz + cz*cz - sy*sy*sz*sz
M₁₂ = -2 * cy * sy * s
M₁₃ =  q_y
M₂₀ = 2 * cy * cz * sy
M₂₁ = -2 * cy * sy * sz
M₂₂ = sy*sy - cy*cy
M₂₃ =  q_z

Thus, 

    | M₀₀  M₀₁  M₀₂  (-q_x*M₀₀ - q_y*M₀₁ - q_z*M₀₂ + q_x) |
R = | M₁₀  M₁₁  M₁₂  (-q_x*M₁₀ - q_y*M₁₁ - q_z*M₁₂ + q_y) |
    | M₂₀  M₂₁  M₂₂  (-q_x*M₂₀ - q_y*M₂₁ - q_z*M₂₂ + q_z) |
    | 0    0    0                        1                |