Physically diffuse reflections are reflections that scatter the incoming light into a multitude of random directions (See: https://en.wikipedia.org/wiki/Diffuse_reflection). However, calculating all of that would require a huge amount of computational resources. It would also result in an exponential growth for the computational complexity as soon as a ray reflects of more than one diffusive surface one after the other. To avoid the exponential growth we only follow a single random ray instead. Approximating the "rest" of the rays will then be dealt with by casting a fix number of rays from the same starting position (Samples per pixel).
Assuming we know the intersection point P_surface as well as the surfaces normal N. We can generate a random ray by using a random angle r1 and a random distance r2 (normalized to be between 0.0 and 1.0). As shown below, we first create an orthonormal coordinate frame (w, u, v) based on the normal N.
Note that u is generated by the cross product of w with either the X or the Y-axis of the default coordinate frame. The decision which is used is made based on the x component of w. If the x component indicates that w is parallel to the Y-Z-plane (|x| < 0.1), the X-Axis is used for the cross product. Otherwise, the Y-Axis is used.
With w, u, and v in place we can now turn our random angle r1 and our random distance r2 into a random vector using this formula: