Generating random directions efficiently is a major issue, especially in GPU implementations of anything like ray tracing. In my research about the matter I've been using an implementation described here. The solution refers to single values but can be easily applied to vectors as well. In my applications I generate random vectors using two methods:
- by repeating calls to PRNG function several times,
- using a function that produces a normally distributed vectors according to given slope value.
In both cases the base functions I use are as follows:
inline __device__ unsigned _TausStep(unsigned &z, const int S1, const int S2, const int S3, const unsigned M)
{
unsigned b = (((z << S1) ^ z) >> S2);
return z = (((z & M) << S3) ^ b);
}
inline __device__ unsigned _LCGStep(unsigned &z, const unsigned A, const unsigned C)
{
return z = (A * z + C);
}
/* Hybrid Tausworthe PRNG by NVIDIA.
https://developer.nvidia.com/gpugems/GPUGems3/gpugems3_ch37.html
*/
inline __device__ float hybtausf(uint4 &z)
{
return 2.3283064365387e-10 * (
_TausStep(z.x, 13, 19, 12, 4294967294UL) ^
_TausStep(z.y, 2, 25, 4, 4294967288UL) ^
_TausStep(z.z, 3, 11, 17, 4294967280UL) ^
_LCGStep(z.w, 1664525, 1013904223UL)
);
}
The hybtausf
function is what I use as a PRNG function. It is used both in the multiple calls random vectors generation, as in the following function, that implements the second method:
// Pseudo-random, normally distributed ray direction generator.
__device__ float3 prrg(float3 dir, float n, float4 &seed)
{
uint4 z = make_uint4(seed.x, seed.y, seed.z, seed.w);
float u0 = hybtausf(z);
float u1 = hybtausf(z);
float3 w = normalize(dir);
float3 t = w;
w = fabs(w);
seed = make_float4(z);
if (w.x < w.y) {
if (w.x < w.z) t.x = 1.0f;
else t.z = 1.0f;
} else if (w.y < w.z) t.y = 1.0f;
else t.z = 1.0f;
w = normalize(dir);
float3 u = normalize(cross(t, w));
float3 v = cross(w, u);
float theta = acos(pow(1.0f - u0, 1.0f / (1.0f + n)));
float phi = 2.0f * PI * u1;
float3 a = make_float3(cos(phi)*sin(theta), sin(phi)*sin(theta), cos(theta));
float3 R[3] = { make_float3(u.x, v.x, w.x), make_float3(u.y, v.y, w.y), make_float3(u.z, v.z, w.z) };
return mult(a, R);
}
The n
parameter controls the mentioned "slope" of the distribution function, while dir
is a "reference direction" around which a new direction will be generated. This method allows to generate the following overall outcome images of my implementation of path tracing (notice the reflection on the sphere):

The image is a bit grainy but that is because I haven't yet applied proper importance sampling in the particular case.
I hope that helps.