When writing illumination related shaders one often needs some directions on the hemisphere oriented at the surface normal. The usual approach is to pre-compute these directions and store them either in a static/uniform array or to create a lookup texture containing the directions. On this page I investigate how reasonably well distributed directions can be quickly computed directly inside the shader using the Hammersley point set in 2d. Of course the point set is not restricted to sampling directions and it could be used for example also for shadow map filtering, screen space ambient occlusion or any other application that needs a 2d sampling pattern.

The Hammersley Point Set in 2d is one of the simplest low discrepancy sequences and has many possible applications in the context of quasi Monte Carlo integration in computer graphics. A general overview of various different quasi Monte Carlo sequences can be found in the book by Niederreiter [niederreiter92]. An explicit use of the Hammersley point set for sampling the sphere can be found for example in [wong97].

Below I first give a short introduction on how to construct the Hammersley point set mathematically and then show how these points look in 2d and how they look when mapped to the hemisphere. Afterwards I show you an implementation that can be directly plugged into a GLSL shader (or easily translated to any other language).

We consider a point set $\mathcal{P_N} = \left\{\ \mathbf{x}_1, ..., \mathbf{x}_N \right\}$ with the number of points $N \geq 1$ on the two-dimensional unit-square $[0, 1)^2$. The Hammersley Point Set $\mathcal{H_N}$ in 2d is now defined as $$\mathcal{H_N} = \left\{\mathbf{x}_i = \binom{i / N}{\Phi_2(i)}, \quad \text{for} \quad i = 0, ..., N - 1 \right\} \qquad\quad (1)$$ where $\Phi_2(i)$ is the Van der Corput sequence.

The idea of the Van der Corput sequence is to mirror the binary representation of $i$ at the decimal point to get a number in the interval $[0, 1)$ The mathematical definition of the Van der Corput sequence $\Phi_2(i)$ is given as: $$ \Phi_2(i) = \frac{a_0}{2} + \frac{a_1}{2^2} + ... + \frac{a_r}{2^{r+1}} \qquad\quad (2)$$ where $a_0 a_1 ... a_r$ are the individual digits of $i$ in binary representation (i.e. $i = a_0 + a_1 2 + a_2 2^2 + a_3 2^3 + ... + a_r 2^r$).

To create now directions on the hemisphere from $\mathcal{H_N}$ any mapping from $[0, 1)$ to the hemisphere could be used but one has to be careful that the chosen mapping fulfills the requirements of the application. When performing integration in the light transport context one usually needs a uniform or cosinus weighted distribution on the sphere.

Let $\mathbf{x}_i = \binom{u}{v} \in \mathcal{H_N}$ be a point from our chosen Hammersley point set. From the two coordinates $u$ and $v$ we now can compute the following mapping.

To implement the Hammersley point set we only need an efficent way to implement the Van der Corput radical inverse $\Phi_2(i)$. Since it is in base 2 we can use some basic bit operations to achieve this. The brilliant book Hacker's Delight [warren01] provides us a a simple way to reverse the bits in a given 32bit integer. Using this, the following code then implements $\Phi_2(i)$ in GLSL:

 float radicalInverse_VdC(uint bits) {
     bits = (bits << 16u) | (bits >> 16u);
     bits = ((bits & 0x55555555u) << 1u) | ((bits & 0xAAAAAAAAu) >> 1u);
     bits = ((bits & 0x33333333u) << 2u) | ((bits & 0xCCCCCCCCu) >> 2u);
     bits = ((bits & 0x0F0F0F0Fu) << 4u) | ((bits & 0xF0F0F0F0u) >> 4u);
     bits = ((bits & 0x00FF00FFu) << 8u) | ((bits & 0xFF00FF00u) >> 8u);
     return float(bits) * 2.3283064365386963e-10; // / 0x100000000
 }

 vec2 hammersley2d(uint i, uint N) {
     return vec2(float(i)/float(N), radicalInverse_VdC(i));
 }

One small optimization would be to replace the division by $N$ in the code above with a multiplication of the inverse since $N$ needs to be fixed anyways.

For the sake of completeness here is also the code to create a uniform or cosinus distributed direction (z-up) from the hammersley point $\mathbf{x}_i = \binom{u}{v}$. Note that the $1 - u$ in the code is necessary to map the first point in the sequence to the up direction instead of the tangent plane.

 const float PI = 3.14159265358979;

 vec3 hemisphereSample_uniform(float u, float v) {
     float phi = v * 2.0 * PI;
     float cosTheta = 1.0 - u;
     float sinTheta = sqrt(1.0 - cosTheta * cosTheta);
     return vec3(cos(phi) * sinTheta, sin(phi) * sinTheta, cosTheta);
 }
    
 vec3 hemisphereSample_cos(float u, float v) {
     float phi = v * 2.0 * PI;
     float cosTheta = sqrt(1.0 - u);
     float sinTheta = sqrt(1.0 - cosTheta * cosTheta);
     return vec3(cos(phi) * sinTheta, sin(phi) * sinTheta, cosTheta);
 }

Here I collect some further notes on how these points could be used and some other things that are good to know.

• 2012-11-07: Initial version of this document

[niederreiter92]  Random number generation and quasi-Monte Carlo methods, Harald Niederreiter , 1992, ISBN:0-89871-295-5

[wong97]   Sampling with Hammersley and Halton Points, Tien-Tsin Wong and Wai-Shing Luk and Pheng-Ann Heng, 1997, Journal of Graphics Tools , vol. 2, no. 2, 1997, pp 9-24.

[warren01]  Hacker's Delight, Henry S. Warren, 2001, ISBN:0201914654

[keller01]   Interleaved Sampling, Alexander Keller and Wolfgang Heidrich, 2001, Proceedings of the Eurographics Workshop on Rendering

[segovia06]   Non-interleaved Deferred Shading of Interleaved Sample Patterns, Benjamin Segovia and Jean-Claude Iehl and Bernard P�roche, 2006, Eurographics/SIGGRAPH Workshop on Graphics Hardware '06