Scrambled Halton

The Halton sequence, which is one of my favourite algorithms ever, can be used for efficient stratified multi-dimensional sampling. Some references:

Van der Corput sequence
Halton sequence

It is possible to do stratified sampling of hyper-points in the s-dimensional unit hyper-cube by picking one consecutive dimension of the Halton series for each component. A convenient way to do so is to use the first s prime numbers as the basis for each Halton sequence.

It is well-known, however, that while this approach works great for low dimensions, high dimensions often exhibit a great degree of undesired correlation. The following image displays a grid where each cell combines two components of the 32-dimensional Halton hyper-cube.

Raw 32-dimensional Halton hyper-cube

One can easily spot how some pairs exhibit an obvious pattern, while others fill their corresponding 2D area very densely. This happens more aggressively for higher dimensions (i.e., to the right/bottom in the image) and for pairs formed with close components (i.e., near the diagonal in the image). Click to enlarge!

A successful approach to dissolve this problem without losing the good properties of stratification is to do “random digit scrambling” (a.k.a. rds). During the construction of a Halton number, digits in the range [0..base[ are combined. Given a Pseudo-Random Permutation of length=base, all that one must do is use PRP(digit) instead of digit directly. This somewhat shuffles Halton pairs in rows and columns in a strong way so the correlation disappears. However, since the PRP is a bijection, the good properties of stratification are generally preserved.

How to build a strong and efficient randomised PRP of an arbitrary length is an interesting subject which details I won’t get into here.

Here’s the scrambling strategy in action:

Scrambled 32-dimensional Halton hyper-cube

Now all the blocks in the grid look equally dense. Magic!

As long as one picks good PRPs, it is possible to generate any number of different samplings, all with the same good properties.

Uniform vs. stratified

This is a classic subject in numerical (Monte Carlo) integration.

Uniform 2D distribution vs. Halton series for the first 2 dimensions

To the left: 32768 points in a 512×512 image using a uniform random number generator (Mersenne Twister). To the right, the first 32768 pairs in the Halton series, using dimensions #0 and #1. Click to enlarge!

Hosek & Wilkie sky model

I have spent the past hours working on the sky model in Arion. This time, I am taking the Hosek & Wilkie model for a spin, as an alternative to our good old implementation of the Preetham sky model:

Preetham et al. sky model

Hosek & Wilkie sky model

Some procedural action

I am (fortunately!) always very busy, but lately it’s been a long chain of little-sleep, endless work sessions. I am working on new technology for Arion 3 and I couldn’t be more excited.

Here’s some image where all the materials are strictly procedural done by Erwann Loison in Arion/ArionFX. Not a single image texture was used:

Procedural action by Erwann Loison

Here is his original post on his blog.

Stuff like this, and much, much more, will be news in the Arion world soon (when it’s done).

Happy New 2015!

Sobel operator

The Sobel operator is a simple way to approximate the gradient of the intensity in an image. This, in visual terms, can be used for edge detection. The purpose of edge detection is to significantly reduce the amount of data in the image, while preserving the structural properties to be used for further image processing.

In practice, the Sobel operator is simply a pair of 3×3 (separable) convolution kernels. One highlights the horizontal gradient/edges, and the other one highlights the vertical gradient/edges.

S_x = \left[ \begin{array}{ccc} -1 & \phantom{+} 0 & +1 \\ -2 & \phantom{+} 0 & +2 \\ -1 & \phantom{+} 0 & +1 \end{array} \right]

S_y = \left[ \begin{array}{ccc} +1 & +2 & +1 \\ \phantom{+} 0 & \phantom{+} 0 & \phantom{+} 0 \\ -1 & -2 & -1 \end{array} \right]

In non-formal terms, and under certain theoretical assumptions, this is conceptually equivalent to computing the partial derivatives of the image with respect to x an y.

For the examples below, I am using the same image featured by Wikipedia in the Sobel operator page:

Sobel operator

This grid presents:

1- The input image (luminance).
2- The absolute magnitude of the result of the Sobel filter.
3- The result of the Sobel (x) filter.
4- The result of the Sobel (y) filter.
5- Same as (2), but in faux color.
6- The gradient vectors, normalized and displayed as an RGB-encoded normal map.

The 3-tap Sobel convolution kernels have a 1px radius, so they have a very limited edge detection ‘range’. This makes the filter quite shaky as soon as the image presents fine details or noise. For this reason, one may want to pre-pass the input image with a subtle Gaussian blur.

This has the effect of diminishing edges in general (as expected), but the resulting gradient images are equivalent, yet much cleaner.

Sobel operator (Gaussian with sigma=2)

The Sobel operator is one of the most fundamental building blocks in feature detection in the field of Computer Vision.

Note that the Sobel operator does not characterize edges or detects features in any way. It simply produces a filtered image where pixels that most likely belong to an area of high gradient (such as an edge) are highlighted.

Bonus remark: Sobel filtering is very similar to what a render engine such as Arion does to transform a height map (bump map) into a normal map.