Handmade Hero » Forums » Code » SampleHemisphere is not uniform
hititduck
1 posts
#13800 SampleHemisphere is not uniform
5 months ago

The current SampleHemisphere function does not generate points distributed uniformly on the hemisphere:

Here's an efficient way to generate uniform samples on the unit sphere attributed to Marsaglia (1972).

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
internal v3
SampleSphere(random_series *Series)
{
    v3 Result = {};
    f32 X1, X2, SquaresSum, SquareRootTerm;

    do
    {
        X1 = RandomBilateral(Series);
        X2 = RandomBilateral(Series);
        SquaresSum = Square(X1) + Square(X2);
    } while (SquaresSum >= 1.0f);

    SquareRootTerm = 2*SquareRoot(1 - SquaresSum);
    Result.x = X1*SquareRootTerm;
    Result.y = X2*SquareRootTerm;
    Result.z = 1 - 2*SquaresSum;

    return(Result);
}


You can flip the sample's position according to its dot product with a normal as before.
Pseudonym
Andrew Bromage
184 posts
1 project

Research engineer, resident maths nerd (Erdős number 3).

#13820 SampleHemisphere is not uniform
5 months ago

Here's a method which involves no loops, due to Archimedes (225 BCE). The idea is to randomly sample the enclosed cylinder then map that onto the sphere or hemisphere.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
internal v3
SampleHemisphere(random_series *Series)
{
    v3 Result = {};
    f32 z = RandomBilateral(Series);
    f32 r = SquareRoot(1 - z*z);
    f32 theta = RandomBilateral(Series) * 2 * PI;
    Result.x = r * Cos(theta);
    Result.y = r * Sin(theta);
    Result.z = z;
    return (Result);
}


However, you probably don't actually want to sample the hemisphere uniformly; for light transport, you probably want to include the geometric cosine factor. The easiest way to do this is to randomly sample the unit disk, then map that onto a hemisphere sitting above it.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
internal v3
CosineSampleHemisphere(random_series *Series)
{
    v3 Result = {};
    f32 r2 = RandomBilateral(Series);
    f32 r = SquareRoot(r2);
    f32 theta = RandomBilateral(Series) * 2 * PI;
    f32 z = SquareRoot(1 - r2);
    Result.x = r * Cos(theta);
    Result.y = r * Sin(theta);
    Result.z = z;
    return (Result);
}

sub f{($f)[email protected]_;print"$f(q{$f});";}f(q{sub f{($f)[email protected]_;print"$f(q{$f});";}f});