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.

The Error function (erf)

Here is the 1D Gaussian function:

G( x , \sigma ) = \frac{ 1 }{ \sqrt{ 2 \pi } \sigma } e^{ \frac{ - x^2 }{ 2 \sigma ^2 } }

Put in short, the Error function is the integral of the Gaussian function from 0 to a certain point x:

\text{erf}( x ) = \frac{ 2 }{ \pi } \int_0^x \! e^{ -t^2 } \, \mathrm{ d }t

At least, that is the way the formula is presented by Wikipedia and Wolfram|Alpha. But as soon as you try to work with it you find out that in order to really match the above Gaussian function, normalization and axis scaling must be taken care of:

\text{erf}( x , \sigma ) = \int_0^{ \frac{ x }{ \sqrt{ 2 } \sigma } } \! e^{ -t^2 } \, \mathrm{ d }t

The plots below display G(x,sigma) (bell-shaped) in blue, and erf(x,sigma) (S-shaped) in yellow.

A very typical use for G(x,sigma) that I’ve been talking about on this blog lately, is to build convolution kernels for Gaussian blur. An image convolution kernel, for example, is a pixel-centered discretization of a certain underlying function (a Gaussian, in this case). Said discretization splits the x axis in uniform unit-length bins (a.k.a. taps, or intervals) centered at x=0.

For each bin, this is pretty much the definition of the integral of G(x,sigma) along the bin. That is, the increment of erf(x,sigma) between both end-points of the bin.

Discretized 1D Gaussian (sigma=0.5)

Discretized 1D Gaussian (sigma=1.0)

Discretized 1D Gaussian (sigma=1.5)

Doing it wrong:

Most implementations of Gaussian convolution kernels simply evaluate G(x,sigma) at the mid-point of each bin. This is a decent approximation that is also trivial to implement:

\text{bin}( a , b , \sigma ) = G( \frac{ a + b }{ 2 } , \sigma )

This value is represented in blue in the above plots.

Implementing erf:

While G(x,sigma) has a trivial explicit formulation, erf is an archetypical non-elementary function.

Some development environments provide an erf implementation. But most (the ones I use) don’t. There are some smart example implementations out there if you Google a bit. Usually they are either numerical approximations, or piece-wise interpolations.

Assuming that one has an erf implementation at his disposal, the correct value for the bin [a..b] is:

\text{bin'}( a , b , \sigma ) = \text{erf}( b , \sigma ) - \text{erf}( a , \sigma )

This value is represented in red in the above plots.

A poor-man approximation:

If you are feeling lazy, an alternative approximation is to super-sample G(x,sigma) along the bin. The amount of samples can be chosen to be proportional to how small sigma is.

This method is easy to implement. However, it is far less elegant and also less efficient as more evaluation calls are required.

Conclusion:

When discretizing a Gaussian function, as sigma becomes small (close to or smaller than (b-a)), the approximation bin(a,b) (in blue above) deviates significantly from the correct bin'(a,b) (in red above).

So if you are going to Gaussian-blur at radius 1px or smaller, and precision is important, then it is necessary to use erf(x,sigma) or at least super-sample G(x,sigma).

Downsampling and Gaussian blur

I talked about several strategies to optimize convolutions in some of my previous posts. I still got to talk about how to approximate a Gaussian blur using a multi-step Box blur in a future post. However, there is yet another good technique to optimize a Gaussian blur that may come handy in some cases.

This post is inspired by a need that I had some days ago: Say that you need to do a 3D Gaussian blur on a potentially humongous 3D data buffer. Working with downsampled data sounds ideal in terms of storage and performance. So that’s what I am going to talk about here:

What happens when downsampling is mixed with a Gaussian blur?

The idea:

Here’s the 0-centered and un-normalized 1D Gaussian function:

G( x , \sigma ) = e^\frac{ - x^2 }{ 2 \sigma ^2 }

The sigma parameter in the Gaussian function stretches the bell shape along the x axis. So it is quite straightforward to understand that if one downsamples the input dataset by a scale factor k, then applies a (smaller) Gaussian where sigma’=s/k, and finally upscales the result by the same scale factor k, the result will approximate a true Gaussian on the original dataset where sigma=s.

In cleaner terms: if one has an input dataset (e.g., an image) I and wants to have it blurred by a Gaussian where sigma=s:

1- I’<=I downsampled by a certain scale factor k.
2- I”<=I' blurred by a small Gaussian where s’=s/k.
3- I”’<=I'' upscaled by a scale factor k.

How good is this approximation?

The Sampling Theorem states that sampling a signal at (at least) twice its smallest wavelength is enough. Which means that downsampling cuts frequencies above the Nyquist limit (half the sampling rate). In other words: Downsampling means less data to process, but at the expense of introducing an error.

Fortunately, a Gaussian blur is a form of low-pass frequency filter. This means that blurring is quite tolerant to alterations in the high part of the frequency spectrum.

Visual evaluation:

In the examples below I am downsampling with a simple pixel average, and I am upscaling with a simple bilinear filter. The 2×2 grids below compare:

1- Top-left – The original image I.
2- Top-rightI downsampled and upscaled by k (note the blocky bilinear filter look).
3- Bottom-left – The resulting image I”’.
4- Bottom-rightI blurred by a true Gaussian where sigma=s.

In these examples, I chose k=sigma for simplicity. This means that the small Gaussian uses sigma’=1.

Gaussian blur where sigma=4

Gaussian blur where sigma=16

Gaussian blur where sigma=64

Conclusion:

As shown, the approximation (bottom-left vs. bottom-right) is pretty good.

The gain in speed depends on multiple implementation factors. However, as I explained above, this post was inspired by a need to cope with a cubic memory storage problem when doing Gaussian blurs on a 3D buffer. Working with a heavily downsampled buffer clearly helps in that sense. And it is needless to say that decreasing the amount of data to process by k^3 also brings a dramatic speed boost, making it possible to use tiny separable convolutions along the 3 (downsampled) axes.

Note that one might pick any downsampling scale factor k>=1. The higher the value of k, the higher the approximation error and the smaller and faster the convolution.

The choice k=sigma offers a good trade-off between approximation error and speed gain as shown above.