Hammersley Points on the Hemisphere


Holger Dammertz

A collection of notes on how to use the Hammersley point set in 2d for fast and practical generation of hemisphere directions in a shader program. If you find any mistakes or have comments do not hesitate to contact me or leave a comment at my blog.

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$).

2.1. Radical Inverse Examples

This is best understood by an example. In the following table I give $i$ in decimal and binary representation and the associated Van der Corput radical inverse (as defined in Equation 2) by mirroring the binary representation at the decimal point. For easier understanding I also include the decimal point in the binary representation of $i$ even though you usually omit this part when considering integral binary numbers.
$i$ decimal $i$ binary $\Phi_2(i)$ binary $\Phi_2(i)$ decimal
0 0000.0 0.0000 0.0
1 0001.0 0.1000 0.5
2 0010.0 0.0100 0.25
3 0011.0 0.1100 0.75
4 0100.0 0.0010 0.125
5 0101.0 0.1010 0.625
11 1011.0 0.1101 0.8125
This table shows the first few numbers in the Van der Corput sequence. It is just the binary representation of $i$ mirrored at the decimal point.

This simplicity of just mirroring the binary representation makes the Hammersley point set an ideal candidate for a quick and easy to implement point set inside a shader. Before I discuss the code I want to show you how the point set actually looks in 2d.

2.2. Image of the Hammersley Point Set

In this interactive image you can change the number $N$ of points and get the resulting point set as defined by Equation 1 visualized in the unit square $[0,1)^2$. As you can see the point $\binom{0}{0}$ is always part of this sequence. You can also see that the points are quite well distributed for any number $N$ of points.

# Points
No Canvas Support

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.

Uniform Mapping

$$ \begin{array}{lcl} \theta & = & cos^{-1}(1 - u) \\ \phi & = & 2 \pi v \end{array} $$

Cosinus Mapping

$$ \begin{array}{lcl} \theta & = & cos^{-1}(\sqrt{(1 - u)}) \\ \phi & = & 2 \pi v \end{array} $$

# Points Mapping
No Canvas Support

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
The $i$th point $\mathbf{x}_i$ is then computed by
 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);
Note that the returned direction is of course in the default coordinate frame with z-up (you can call it also tangent space if you like). So the last step that probably needs to be performed in an application is transforming the returned direction into world space by using a ortho-normal basis matrix constructed from your local normal vector.

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

5.1. How good are these points?

Quality is of course application dependent but the first important thing to note is that the Hammersley point set in 2d is not very nicely distributed on the hemisphere for a small number of points. This can be easily verified above in the section <Generating Points on the Hemisphere>. This problem lessens when many points are needed as for example when interleaved sampling is performed. While there are many other quasi-random sequences that could be used, so far I haven't found any that can be implemented as efficiently as the Hammersley point set (except for rank-1 lattices but they can't be used easily for varying number of points).

Another aspect of this sampling pattern is that the number $N$ of points needs to be known a-priori before starting to compute any point. This is almost never a problem in real time rendering applications, but as soon as some form of progressive rendering is performed using a quasi-monte carlo sequence (instead of a point set) is probably a better choice.

5.2. Interleaved Sampling

Interleaved sampling [keller01] is an often used optimization [segovia06] for real time rendering where the computation of an integral is spread over several pixel and later combined in an additional filtering step.

Of course any collection of points can be used for this but the easy construction and usage makes it espacially simple to split the point set over several fragments and/or hemispheres.

  • 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]  pdf 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]  pdf Interleaved Sampling, Alexander Keller and Wolfgang Heidrich, 2001, Proceedings of the Eurographics Workshop on Rendering

[segovia06]  pdf 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