I messed up this formula in my course this week. I promised to fix it here. Here it is!
Suppose you want to emit ray from a Lambertian (ideal diffuse) surface. This has the probability of a direction proportional to cosine of the angle with the normal. If we use a spherical coordinate system with (theta,phi) where phi is the angle "around" the z-axis and theta is the angle "down from" the z axis, we can see from symmetry that all phi are equally likely so
phi = 2*PI*random(0,1)
Now theta you can't just treat as a 1D random variable with density proportional to cosine of theta. The "differential measure" on the sphere is sin(theta)*dtheta*dphi. The sine of theta goes along with the theta term so the pdf of theta as a 1d density function is proportional to sin(theta)cos(theta)
p(t) = k*sin(t)*cos(t) // t is theta
We know INT p(t) dt = 1, with t ranging 0 to PI/2 (hemisphere). We will have to compute the cumulative probability density function and can normalize then. The key bit is this: if we call r = random(0,1), we can get a Theta by solving for:
r = INT_0^Theta p(t) dt
The basic idea behind that is that if r = 0.7 for example, we want 70% of the "mass" of the density function to be below the Theta we choose.
r = INT_0^Theta k*sin(t)*cos(t)
r = 0.5*k*sin^2(Theta)
given Sin goes 0 to 1 the normalization falls out easily:
sin^2(Theta) = r
In spherical coordinates we have z = cos(Theta), so a unit vector with cosine density has the Cartesian coordunates:
x = cos(2*PI*r1)*sqrt(r)
y = sin(2*PI*r2)*sqrt(r)
z = sqrt(1-r)