Graphics Programming and Theory - Articles

Dynamic 2D Soft Shadows

Fri, 19 Apr 2013 12:42:15 +0000

The aim of this document is to describe an accurate method of generating soft shadows in a 2D environment. The shadows are created on the fly via geometry as opposed to traditional 2D shadow methods, which typically use shadow sprites or other image based methods. This method takes advantage of several features available in the core OpenGL API, but is by no means restricted to this platform. Direct3D would also be a suitable rendering API, and the concepts and reasoning behind the various rendering stages should hopefully be clear enough that a reader can make the conversion without too much hassle.

Note:
This article was originally published to GameDev.net back in 2004. It was revised by the original author in 2008 and published in the book Advanced Game Programming: A GameDev.net Collection, which is one of 4 books collecting both popular GameDev.net articls and new original content in print format.

Overview

We will start by defining a few terms that we will use frequently, and a brief explanation of the phenomena that we are attempting to reproduce on our digital canvas.

Image 1: Overview of terms

Light source

An obvious place to start – in this implementation we will discuss a point light source, although extending the method to include directional lights as well would be easily done, as is adding ambient lighting into the system. We use a point light source with a user-defined radius to generate the soft shadows accurately.

Shadow caster

A shadow caster is any object that blocks the light emitted from the source. In this article we present implementation details for using convex hulls as shadow casters. Convex hulls have several useful properties, and provide a flexible primitive from which to construct more complex objects. Details of the hulls are discussed in just a bit.

Light range

In reality light intensity over a distance is subject to the inverse square relationship, and so can never really reach zero. In games however linear light fall off often looks as good or better depending on the circumstances. The image above a linear fall off in the intensity is used, dropping to zero at the edge of the light range.

Umbra

The umbra region of a shadow is the area completely hidden from the light source, and as such is a single colour (the image above shows the umbra region in black since there is no other light source to illuminate this region).

Penumbra

The penumbra region of a shadow is what creates the soft edges. This is cast in any area that is partially hidden from the light but neither in full view or totally hidden. The size and shape of the penumbra region is related to the lights position and physical size (not its range).

Core Classes

First we'll have a look at a couple of classes that are at the core of the system – the Light and ConvexHull classes.

Light: The light class is fairly self-explanatory, holding all the information needed to represent a light source in our world.

Contains:

Position and depth. Fairly obvious, these are the location in the world. Although the system is 2d, we still use a depth for correctly defining which objects to draw in front of which others. Since we'll be using 3D hardware to get nice fast rendering we'll take advantage of the depth buffer for this.
Physical size and range. Both stored as a simple radial distance, these control how the light influences its surroundings.
Colour and intensity. Lights have a colour value stored in the standard RGB form, and an intensity value which is the intensity at the centre of the light.

ConvexHull: The convex hull is our primitive shape from which we will construct our world. By using these primitives we are able to construct more complex geometry.

Contains:

List of points. A simple list is maintained of all the points that make up the edges of the hull. This is calculated from a collection of points and the gift-wrapping algorithm is used to discard unneeded points. The gift-wrapping method is useful since the output geometry typically has a low number of edges. You may want to look into the QuickHull method as an alternative.
Depth. As for the light, a single depth value is used for proper display of overlapping objects.
Shadow depth offset. The importance of this is described later.
Centre position. The centre of the hull is approximated by averaging all the points on the edge of the hull. While not an exact geometric centre it is close enough for our needs.
Vertex data. Other data associated with the vertex positions. Currently only per-vertex colours but texture cords could be added without requiring any major changes.

Rendering Overview

The basic rendering process for a single frame looks like:

Clear screen, initialise camera matrix
Fill z buffer with all visible objects.
For every light:
1. Clear alpha buffer
2. Load alpha buffer with light intensity
3. Mask away shadow regions with shadow geometry
4. Render geometry with full detail (colours, textures etc.) modulated by the light intensity.

The essential point from the above is that a rendering pass is performed for every visible light, during which the alpha buffer is used to accumulate the lights intensity. Once the final intensity values for the light have been created in the alpha buffer, we render all the geometry modulated by the values in the alpha buffer.

Simple Light Attenuation

First we'll set up the foundation for the lighting – converting the above pseudo code into actual code but without the shadow generation for now.

public void render(Scene scene, GLDrawable canvas)
{
GL gl = canvas.getGL();
gl.glDepthMask(true);
gl.glClearDepth(1f);
gl.glClearColor(0.0f, 0.0f, 0.0f, 0.0f);
gl.glClear(GL.GL_COLOR_BUFFER_BIT |
GL.GL_DEPTH_BUFFER_BIT |
GL.GL_STENCIL_BUFFER_BIT);
gl.glMatrixMode(GL.GL_PROJECTION);
gl.glLoadIdentity();
gl.glMatrixMode(GL.GL_MODELVIEW);
gl.glLoadIdentity();
gl.glMatrixMode(GL.GL_TEXTURE);
gl.glLoadIdentity();
gl.glDisable(GL.GL_CULL_FACE);
findVisibleLights(scene);
Camera activeCamera = scene.getActiveCamera();
activeCamera.preRender(canvas);
{
// First we need to fill the z-buffer
findVisibleObjects(scene, null);
fillZBuffer(canvas);
// For every light
for (int lightIndex=0; lightIndex
Note that code here is written in Java, using the Jogl set of bindings to OpenGL. For C++ programmers you simply have to remember that primitives such as int, float, boolean etc. are always passed by value, and objects are always passed by reference. OpenGL commands and enumerations are scoped to a GL object, which leads to the slightly extended notation from the straight C style.

First we reset the GL state ready for the next frame, collect all the lights that we will need to render this frame and retrieve the currently active camera from the scene. Camera.preRender() and .postRender() are used to set the modelview and projection matrices to that needed for the view position.

Once this initialisation is complete we need to fill the z-buffer for the whole scene. Although not discussed here, this would be the perfect place to take advantage of your favourite type of spatial tree. A quad-tree or AABB-tree would make a good choice for inclusion within the scene, and would be used for all testing of objects against the view frustum. To fill the depth buffer we simply enable z-buffer reading and writing, but with colour writing disabled to leave the colour buffer untouched. This creates a perfect depth buffer for us to use and stops later stages blending pixels hidden from view. It is worth noting that by enabling colour writing an ambient lighting pass can be added here to do both jobs at the same time. From this point onwards we can disable depth writing as it no longer needs to be updated.

Now we perform a rendering pass for every light.

First the alpha buffer is cleared in preparation for its use. This is simply a full screen quad drawn without blending, depth testing or colour writing to reset the alpha channel in the framebuffer to 0f. Since we don't want to disturb the current camera matrices that have been set up, we create this quad by using the current camera position to determine the quads coordinates.

Next we need to load the lights intensity into the alpha buffer. This does not need any blending, but depth testing is enabled this time to allow lights to be restricted to illuminating only the objects beneath them. Again colour writing is left disabled since we are not ready to render any visible geometry yet. The following function is used to create the geometry for a single light:

public void renderLightAlpha(float intensity, GLDrawable canvas)
{
assert (intensity > 0f && intensity <= 1f);
GL gl = canvas.getGL();
int numSubdivisions = 32;
gl.glBegin(GL.GL_TRIANGLE_FAN);
{
gl.glColor4f(0f, 0f, 0f, intensity);
gl.glVertex3f(center.x, center.y, depth);
// Set edge colour for rest of shape
gl.glColor4f(0f, 0f, 0f, 0f);
for (float angle=0; angle<=Math.PI*2;
angle+=((Math.PI*2)/numSubdivisions) )
{
gl.glVertex3f( radius*(float)Math.cos(angle) + center.x,
radius*(float)Math.sin(angle) + center.y, depth);
}
gl.glVertex3f(center.x+radius, center.y, depth);
}
gl.glEnd();
}

What happens is we create a triangle fan rooted at the centre position of the light, then loop around in a circle creating additional vertices as we go. The alpha value of the centre point is our light intensity, fading linearly to zero on the edges of the circle. This creates the smooth light fall off seen in the first image. If other methods of light attenuation are needed, they can be generated here. An interesting alternative would be to use an alpha texture instead of vertex colours; a 1D texture could happily represent a non-linear set of light intensities. Other unusual effects could be achieved by animating the texture coordinates over a 2D texture, such as flickering candlelight or a pulsing light source.

So now we have our light intensity values in the alpha buffer, we will skip the generation of shadow hulls for the moment and move on to getting our level geometry up on the screen.

The geometry pass is where we really start to see things coming together, using the results we have carefully composed in the alpha of the framebuffer. First we need to make sure we have depth testing enabled (using less-than-or-equals as before), and then enable and set up our blending equation correctly.

gl.glEnable(GL.GL_BLEND);
gl.glBlendFunc(GL.GL_DST_ALPHA, GL.GL_ONE);

Simple, yes? What we're doing here is multiplying our incoming fragments (from the geometry we're about to draw) by the alpha values already sitting in the framebuffer. This means any alpha values of 1 will now be drawn at full intensity, and values of 0 being unchanged. This is then added to the current framebuffer colour multiplied by one. This addition to the existing colour means we slowly accumulate our results from previous passes. With our blend equation set up, we simply render our geometry as normal, using whatever vertex colours and textures that takes our fancy.

If you take another look at our render() function near the top, you'll see we've almost finished composing our frame. Once we've looped over all the lights we've practically finished, but we'll insert a couple of extra stages. First an emissive or self illumination pass – this is discussed near the end of the article. After this is a simple wireframe rendering with draws object outlines such as seen in the first image.

Image 2: Per pixel lighting with intensities accumulated in the alpha buffer.

Coloured Lighting

What was once seen as 'the next big thing' in the Quake 2 and Unreal era, coloured lighting is pretty much standard by now, and a powerful tool for level designers to add atmosphere to a scene. Now since we've already got our light intensity ready and waiting for our geometry in the alpha buffer, all we need to do is modulate the geometry colour by the current light colour while drawing. That's a whole lot of multiplication if we want to do it ourselves, but on TnL hardware we can get it practically for free with a simple trick. We enable lighting while drawing our geometry; yet define no normals for we have no need of them. Instead we just enable a single light and set its ambient colour to the colour of our current light. The graphics card will calculate the effect of the light colour on our geometry for us and we need barely lift a finger. Note that because we're accumulating light intensities over multiple lights in the framebuffer we get accurate over brightening effects when lights overlap, and multiple coloured lights will merge and produce white illumination of our objects.

Hard-edged Shadow Casting

Now we have our lights correctly illuminating their surroundings we can start thinking about correctly limiting their light to add shadows into the scene. First we will cast hard edged shadows from shadow casters and then extend this to cover soft edged shadows with correct umbra and penumbra. This is done in the function we previously skipped, mergeShadowHulls().

You will remember that at this point in the rendering we have the light intensity stored in the alpha buffer. Now what we will do is create geometry to represent the shadow from each shadow caster, then merge this into the alpha buffer. This is done inside the ConvexHull class.

Finding the boundary points

Our first step is to determine which points our shadow should be cast from. The list of points that make up the ConvexHull is looped though, and each edge is classified in regard to the light position. In pseudo code:
For every edge: Find normal for edge<
Classify edge as front facing or back facing
Determine if either edge points are boundary points or not.
The normal for the edge is found as:

float nx = currentPoint.y - prevPoint.y;
float ny = currentPoint.x - prevPoint.x;

Then a dot product is performed with this vector and the vector to the light. If this is greater than zero, the edge is front facing. Once and edge has been classified, it is compared against the previous edge. If one is front facing and the other back facing, then the shared vertex is a boundary point. As we walk around the edge of the hull (in an anti-clockwise direction) the boundary point from light to shadow is the start shadow point. The boundary from shadow to light is the end shadow point.

Creating the Shadow Geometry

Once we have these positions, we can generate our shadow geometry. Since we are only generating hard edged shadows at the moment, we will be ignoring the physical size of our light source. Image 3 shows a how the shadow geometry is built.

Image 3: Hard-edged shadow generation

As shown in the image, the shadow geometry is a single triangle strip projected outwards from the back facing edges of the shadow caster. We start at the first boundary point (marked with a red cross) and work our way anti-clockwise. The second point is found by finding the vector from the light to the point, and using this to project the point away from the light. A projection scale amount is used to ensure that the edges of the shadow geometry are always off screen. For now we can simply set this to a sufficiently large number, but later it will be advantageous to calculate this every frame depending on how far zoomed in or out the camera is.

We render the shadow geometry with depth testing enabled to properly layer the shadow between various other objects in the world, but with colour writing disabled, only the alpha in the framebuffer is changed. You may remember that the final geometry pass is modulated (multiplied) by the existing alpha values, which means we need to set the shadow to have an alpha value of zero. Because the framebuffer will clamp the values to between one and zero, overlapping shadows will not make an affected area too dark but instead merge correctly.

Image 4: Hard-edged shadows

Notice in image 4 how the shadow from the diamond correctly obscures the rightmost object, and that their shadows are correctly merged where they overlap.

Soft-Edged Shadow Casting

Now we can properly construct hard edged shadows it is time to extend this to cover soft shadows – note that we cannot simply add faded edges to the existing shadow geometry, since this would result in inaccurate penumbra and umbra regions. First we start by defining a physical radius for the light source to generate the correct penumbra regions, then we need to create the penumbra geometry and modify the creation of the umbra region that we used for the hard edged shadows.

Shadow Fins

Each penumbra region will be created by one or more shadow fins via that ConvexHull and ShadowFin classes.

ShadowFin: An object to encompass all or part of a penumbra region.

Contains:
Root position. This is the position from which the fin protrudes from.
Penumbra vector. This is a vector from the root position which lies on the outer edge of the fin (the highest light intensity).
Umbra vector. This vector from the root position lies on the inner edge of the fin (lowest light intensity).
Penumbra and umbra intensities. These are the light intensities of their relative edges for the fin. If the fin makes up an entire penumbra region these are one and zero respectively.
We start at the first boundary point, and create a ShadowFin from this point. The root position becomes the boundary point, and the penumbra and umbra intensities are initially one and zero. The difficult part of the fin – the penumbra and umbra vectors – is done by the getUmbraVector and getPenumbraVector methods within our Light object.

Image 5: Shadow fin generation

If we look at the vector that lies along the outer penumbra edge we can imagine it as the vector from the light though the boundary point (C, the centre vector) displaced by the light radius. So we must find this displacement.

First we note that the displacement is as right angles to the centre vector. So we take C and find this perpendicular vector in the same way we did to find the normals for the hull edges. Now although looking at the diagram we know which way we want this to point, when we're dealing with boundary points and light positions at all sorts of funny angles to each other we may end up with it pointing in the opposite direction to that which we expect. To solve this we find the vector from the centre of the hull to the boundary point (between the two Xs in the image), and take the dot product of this and the perpendicular vector. If this is less than zero, our vector is pointing in the wrong direction and we invert it.

Armed with this new vector we normalise it and the centre vector, then add them together and we've found our crucial outer penumbra vector. Finding the inner vector requires we repeat the process but this time we invert the logic for the dot product test to displace the centre vector in the opposite direction. We now have a fully calculated shadow fin to send to our renderer!

Non-Linear Shading

Although we have all the numbers we need to render our shadow fin, we'll soon hit a major snag – we can't use regular old vertex colours this time to write to the alpha buffer. We need the inner edge of the penumbra to be zero intensity (zero alpha) and our outer edge to be fully bright (alpha of one). While you can probably visualise that easily, getting our graphics card to actually render a triangle with the colours like this just isn't possible. Try it yourself if you're not sure, you'll soon see how it's the root point that causes the problems – it lies on both edges, so needs to be 'both' zero and one at the same time.

The solution to this (and indeed most cases when you need non linear shading) is to abandon vertex colours for the fins and instead use a texture to hold the information. Below is a copy of the texture I used.

Image 6 : Penumbra texture

You can clearly see how the shadow fin will be rooted at the bottom left, and the two edges running vertical and diagonally to the top edge. Since we don't want texels from the right edge bleeding into the left we set the texture wrapping mode to clamp to the edge values (using glTexParameteri and GL_CLAMP_TO_EDGE). The bottom right half of the texture is unused, although if you really wanted to you could pack something else in here just as long as you're careful not to bleed over the edge.

So we load this texture and bind it for use before drawing our shadow fins, and set the vertex colour to white to leave the texture unchanged by it. Other than that rendering the fins is no different from the shadow hull. The only other thing we need to watch out for is how far back we project our points by the umbra/penumbra vectors, as the limited resolution of our penumbra texture will show if these are moved too far away. Ideally they will be projected to just off screen.

Modifying the umbra generation

Now we've got the fins drawn, we can fill in the umbra between them. This is done in almost exactly the same way as with hard shadows, except we must use the fins inner edges to start and finish from instead of directly projecting away from the centre of the light source. As we move across the back of the shadow caster, we perform a weighted average between these two edge vectors to properly fill in the umbra region. When done correctly we see no gaps between the fins and the umbra geometry, giving us one consistent, accurate shadow cast into the alpha buffer.

Making it robust

Self-Intersection

Now we have this far implemented, the shadows will be looking quite nice – when static – however problems will become apparent when moving the light sources around. The most glaring is that of self-intersection.

Image 7: Self intersection of shadow fin

If the light source is too large in relation to the object, or too near, the inner penumbra edge will intersect the hull. First we need to detect this case. We find the vector from the boundary point to the next edge point in the shadow (moving anti clockwise here since we're on the shadow start boundary point). Then we compare the angle between the outer penumbra edge and our newly found hull edge, and the angle between the outer and inner penumbra edges. If the angle to the boundary edge is smaller, then we've got an intersection case we need to fix.

First we snap the current fin to the boundary edge, and calculate the new intensity for the inner edge via the ratio of the angles. Then we create a new shadow fin at the next point on the hull. This has an outer edge set to the same vector and intensity as the existing fins inner edge, while the new fins inner edge is calculated as before. By gluing these two fins together we create a single smooth shadow edge. Technically we should repeat the self-intersection test and keep adding more fins as needed, however I've not found that this is needed in practice.

Eliminating 'Popping'

You will also notice one other problem with this as it stands, the shadow fins will 'pop' along the edges of the hull as a light rotates. This is because we're still using an infinitely small point light to find the boundary points. To solve this we should take the physical radius into account when finding them. A robust way of doing this is to shift the light source position towards the hull by the radius distance before we find our boundary points. With these two fixes in place the fins will be visually stable as either the light or the hull moves (or both!).

Depth Offset

Depending on the style of game and the view used (such as a side scrolling platformer as opposed to a top down shooter) the way light and shadow interacts with the various level objects will be different. What seems sensible for one may appear wrong in another. Most obviously is with objects casting shadows onto objects at the same depth.

Image 8: The effect of shadow offset

The image above shows the same scene with different shadow offsets. Imagine that the scene is a top down viewpoint – the light grey areas are impassable walls surrounding the floor showing a T junction (imagine hard!). Now the image on the right seems slightly out of place – the shadows are being projected on top of the walls, yet these are much higher than the light source – realistically they shouldn't be lit at all but solid black walls aren't very visually appealing. The second shows the shadows being offset and only obscuring the floor and any objects on it.

Now if you were to imagine the same scene as a 2D platformer, you might prefer the left image. Here it seems to make more sense that the objects should shadow those on the same level. This decision is usually very dependent on the geometry and art direction of the level itself, so no hard and fast rules seem to apply. The best option seems to be to experiment and see which looks best.

Adding control over this is a small modification. At the moment the scene on the left is the common case, and by generating shadow volumes that are a close fit to the edge of the shadow caster we've already done all the hard work, all we need to do is store a shadow depth offset in our ConvexHull and apply it to the depth of the shadow geometry. The existing depth testing will reject fragments drawn behind other objects and leave them at the original intensity.

Emmissive / Self Illumination pass

This is a simple addition that can produce some neat effects – and can be seen as a generalisation of the wireframe 'full-bright' pass. After the lights have been drawn, we clear the alpha buffer again as before, but instead of writing light intensities into it we render our scene geometry with their associated emissive surface. This is an alpha texture used to control light intensities as before, and can be used for glowing objects, such as a computer display or a piece of hardware with a bank of LEDs - everything that has its own light source but is too small to require an individual light of its own. Because these are so small, we skip the shadow generation and can do them all in one go. Then we modulate the scene geometry by this alpha intensity as before. Unusual effects are possible with this, such as particles which light up their immediate surroundings, or the bright local illumination provided by a neon sign (with one or two associated lights to provide the lighting at medium and long range).

Scissor Testing

We are extending the shadow geometry until it's off the edge of the screen, but often the area a light affects is much smaller than this. The scissor test (glScissor in OpenGL) allows us to restrict rendering to a rectangle within our window and avoid drawing pixels that have no effect. We just have to project the light's bounds to screen space and set the scissor area before drawing the shadow geometry. This can increase the framerate considerably.

Conclusion

After a lot of work, much maths and a few sneaky tricks and we've finally got a shadow system that's both physically accurate and robust. The hardware requirements are modest indeed – a framebuffer with an alpha component is about all that's required, we don't even need to stray into extensions to get the correct effect. There is a whole heap of possible optimisations and areas for improvement, notably caching the calculated shadow geometry when neither the light nor the caster has changed, and including some sort of spatial tree so that only the shadow casters that are actually needed are used for each light.

Thin Film Interference for Computer Graphics

Fri, 08 Mar 2013 03:59:03 +0000

Note

The theoretical parts of this article require some knowledge of optics and electromagnetism, however the conclusion and final result (a practical implementation of single-layer thin film interference in the context of a BRDF) do not. You may therefore wish to skip the theoretical sections.

Introduction

Wave interference of light has been neglected for a long time in computer graphics, for multiple reasons. Firstly, it is often insignificant and can be cheaply approximated, or even ignored completely. Secondly, it is harder to understand as it requires interpreting light as waves instead of particles (photons). However, interference crops up almost everywhere in daily life, and has recently gained popularity in rendering applications.

Examples of wave interference of light are soap bubbles, gasoline rainbow patterns, lens flares, basically everything that looks cool and/or involves multicolor patterns. For instance, in computer graphics, soap bubbles were in the past approximated with more or less realistic multicolor textures slightly panned with view angle. But it turns out that they are not that computationally difficult to accurately render. We will learn how.

This article will focus on one particular form of interference, namely thin film interference. This occurs when one or more very thin transparent coatings ("films") are placed on top of a material. The films are so thin that when a light wave comes into contact with these film layers, it reflects and refracts multiple times inside the layer system, and interferes with itself in the process.

The goal is to calculate the amount of light reflected off the layer system, and the amount of light transmitted into the internal medium. We will make the assumption that no light is absorbed, which is not required but makes the calculations more approachable as considering absorption of light involves delving deep into Maxwell's equations (the behaviour of electromagnetic waves at interfaces of lossy media is nontrivial). Though in general, each layer is so thin that absorption effects can be neglected most of the time.

We will derive a physical solution for the case where only one film is present (single-layer) and conclude on how to solve the general case with arbitrarily many layers. The single-layer case is sufficient to render most real life occurrences of thin-film interference, however using more layers enables many more advanced effects. The cost of calculating reflection and transmission coefficients is linear in the number of layers.

Derivation

Consider a light wave incident to a thin layer of depth and real refractive index . The external medium has refractive index and the internal medium has refractive index . The incident angle made by the incident light wave and the film's surface normal is , the angle inside the layer is and the refracted angle (inside the internal medium) is denoted . We will also give numbers to each of the three media: medium 0 is the external medium, medium 1 is the layer, and medium 2 is the internal medium. Also, naturally, medium 0 has to have a different refractive index than medium 1, and the same goes for medium 1 and medium 2. Media 0 and 2 can be the same, of course.

First, we know from Snell's Law that the following holds:

Therefore the angles and can be derived from .

Now, we see from the diagram that the only path the light wave can follow is a zigzag pattern as it bounces back and forth between the layer, until it gets transmitted either back into the external medium or into the internal medium. We also note that all the reflected waves (denoted , , ...) and all the transmitted waves are parallel. This is necessary for interference and is a natural consequence of the reciprocal nature of Snell's Law.

And we have assumed that the media involved are non-absorbing, therefore by conservation of energy the reflection and transmission coefficients must sum up to exactly one. It turns out that it is slightly easier to derive the transmission coefficient, so we will do that, but we would get the same thing either way. The reason for this is because the very first reflected wave does not actually penetrate the layer, which means it needs to be handled separately. This does not occur for transmitted waves.

Now let's take a look at what happens to the amplitude of the light wave as it travels through this layer system. First, we need to introduce the Fresnel equations, which let us calculate how much of a light wave's amplitude is reflected and how much of it is transmitted whenever it comes into contact with an interface. These equations should be familiar, although perhaps not in the following form:

These are amplitude reflection/transmission coefficients, for s-polarized and p-polarized light. Indeed, light polarization is important, and in the derivation we will assume the light wave has a given known polarization.

We will now introduce some notation. The following denotes the amplitude reflection coefficient for a light wave going from medium to medium :

Where the correct reflection coefficient is chosen based on the light wave's polarization. Similarly, the amplitude transmission coefficient is:

Because the refractive indices and incident angles for each medium are known and constant, we do not need to specify them.

We are now ready to tackle the problem. Consider the transmitted wave . It's easy to see that since it crosses the layer at two locations, and never reflects anywhere, its amplitude will be:

What about the second transmitted wave ? This one is transmitted once from medium 0 to medium 1, reflects off the medium 1 to medium 2 interface, is reflected again from the medium 1 to medium 0 interface, and is finally transmitted across the medium 1 to medium 2 interface. So its amplitude will be:

We can see there's a pattern here. Every successive transmitted wave will simply reflect two additional times off the top and bottom interface. So, if we denote the amplitude of the th transmitted wave , we have:

We note that even though there are (in theory) infinitely many transmitted waves, their amplitude decreases exponentially, since the Fresnel amplitude reflection coefficient is never quite 1 (except in the case of total internal reflection, where all light is reflected and none is transmitted, of course, if this is the case then the incident wave fully reflects off the layer first chance it gets and so this analysis doesn't apply).

We now have the amplitudes of each transmitted wave. Can we calculate the total amount of transmitted light now? Not quite. These are waves, and you can't just add waves using their amplitudes. We need to consider the phase of each transmitted wave, as these waves might cancel each other out depending on their phase (out of phase waves cancel out, in phase waves amplify each other). The waves also have a frequency, but the frequency depends only on the incident wave's wavelength, which is known and constant, so it can be taken out of the equation.

How do we calculate the phase of each transmitted wave? This is in fact a simple textbook thin film interference problem, and if we denote the phase of the th transmitted wave , the following holds:

Where is the light wave's wavelength and is a constant meant to account for phase changes upon reflection (we will expand on this soon). The important thing is that the phase of every transmitted wave is a multiple of a constant (with respect to the wave index )! That is:

The explanation for this lies in the rather trivial observation that the distance travelled by the light wave inside the layer increases by a constant factor for every consecutive transmitted wave. This is fortunate, as it makes the upcoming calculations very simple. Had the phase depended on in a more complicated way, the problem could have very well been analytically intractable.

We will now explain the meaning of the term. When a wave (any wave, not just electromagnetic light waves) reflects off a medium denser than the one it is in, it will undergo a 180-degree phase change. Because the refractive index is a measure of how dense a medium is, we can use that to calculate this constant. There are two possible reflections here: one at the top interface and one at the bottom interface. We denote:

For the reflection phase change when reflecting off the interface from medium to medium . Therefore, we see that:

Which is constant, as it depends only on the refractive indices of each medium.

At this point we have the amplitude and phase of each transmitted wave. All we have to do is sum them up (as waves), and take the squared magnitude of the resulting complex amplitude to obtain the transmitted intensity. However, because the transmitted waves are in a different medium than the incident wave, we need to take into account the ratio of beam surface area to make sure energy is conserved. That is, we need to multiply by:

This is actually two factors in one. The first, ratio of refractive indices, is there because the transmitted wave won't, in general, have the same speed as the incident wave (for instance, light travels slower in water than in air). So the perceived intensity will not be the same. Remember, intensity is energy per second per squared area, so if the wave is faster the intensity will be higher, so we need to scale the intensity down by a corresponding amount to make sure energy is conserved. The second factor, ratio of cosines, exists because of the change in area of a beam of light as it is refracted. The following diagram illustrates all of this nicely:

It is worth noting that reflected light is treated the same, however because reflected waves remain in the same medium and the reflected angle is the same as the incident angle, both ratios just cancel out.

Now, we have the following expression for the transmitted intensity:

This looks complicated, but it actually isn't. This is because both the phase and the amplitude are dependent on in such a way that:

And we will now use the following two substitutions, just to make the expressions a bit more readable:

We now have a geometric series sum, which we can evaluate as follows:

Simplifying rather elegantly to the following (assuming is real):

And by conservation of energy, we have:

Which concludes the derivation. As a final note, we can calculate the average over all possible phases of this result. If we are correct, then we should get the same result as a geometric optics derivation. The reader can verify that, indeed, we have:

We conclude that the transmission coefficient (intensity of light transmitted across the layer) is and the reflection coefficient (intensity of light reflected off the layer) is .

General Case

The derivation shown above is quite naive, and does not generalize well at all to multiple layers, though it is the simplest method to see what is happening at a low level. If you wish to implement n-layer thin film interference, the method of choice is the Transfer-matrix method, which simplifies the problem down to a series of matrix multiplications and can be derived using powerful electromagnetism techniques.

Implementation

So we now know just how much light is reflected from the layer. How can we implement this in the context of a BRDF? It's quite simple: this reflected term simply replaces the ordinary Fresnel term, accounting for thin film interference effects. This means you can trivially include thin-film interference effects in any BRDF as long as it has a Fresnel term. The function below computes the reflection coefficient for a given wavelength and incident angle.

// cosI is the cosine of the incident angle, that is, cos0 = dot(view angle, normal)
// lambda is the wavelength of the incident light (e.g. lambda = 510 for green)
float ThinFilmReflectance(float cos0, float lambda)
{
    const float thickness; // the thin film thickness
    const float n0, n1, n2; // the refractive indices
    
    // compute the phase change term (constant)
    const float d10 = (n1 > n0) ? 0 : PI;
    const float d12 = (n1 > n2) ? 0 : PI;
    const float delta = d10 + d12;
    
    // now, compute cos1, the cosine of the reflected angle
    float sin1 = pow(n0 / n1, 2) * (1 - pow(cos0, 2));
    if (sin1 > 1) return 1.0f; // total internal reflection
    float cos1 = sqrt(1 - sin1);
    
    // compute cos2, the cosine of the final transmitted angle, i.e. cos(theta_2)
    // we need this angle for the Fresnel terms at the bottom interface
    float sin2 = pow(n0 / n2, 2) * (1 - pow(cos0, 2));
    if (sin2 > 1) return 1.0f; // total internal reflection
    float cos2 = sqrt(1 - sin2);
    
    // get the reflection transmission amplitude Fresnel coefficients
    float alpha_s = rs(n1, n0, cos1, cos0) * rs(n1, n2, cos1, cos2); // rho_10 * rho_12 (s-polarized)
    float alpha_p = rp(n1, n0, cos1, cos0) * rp(n1, n2, cos1, cos2); // rho_10 * rho_12 (p-polarized)
    
    float beta_s = ts(n0, n1, cos0, cos1) * ts(n1, n2, cos1, cos2); // tau_01 * tau_12 (s-polarized)
    float beta_p = tp(n0, n1, cos0, cos1) * tp(n1, n2, cos1, cos2); // tau_01 * tau_12 (p-polarized)
        
    // compute the phase term (phi)
    float phi = (2 * PI / lambda) * (2 * n1 * thickness * cos1) + delta;
        
    // finally, evaluate the transmitted intensity for the two possible polarizations
    float ts = pow(beta_s) / (pow(alpha_s, 2) - 2 * alpha_s * cos(phi) + 1);
    float tp = pow(beta_p) / (pow(alpha_p, 2) - 2 * alpha_p * cos(phi) + 1);
    
    // we need to take into account conservation of energy for transmission
    float beamRatio = (n2 * cos2) / (n0 * cos0);
    
    // calculate the average transmitted intensity (if you know the polarization distribution of your
    // light source, you should specify it here. if you don't, a 50%/50% average is generally used)
    float t = beamRatio * (ts + tp) / 2;
    
    // and finally, derive the reflected intensity
    return 1 - t;
}

We can now sample this function at red, green, and blue wavelengths (650, 510, 475 nanometers, respectively) and substitute the RGB reflectance obtained into the Fresnel term of the BRDF. Or, if you are rendering spectrally, just give the wavelength directly. That's it.

One word on polarization - in general, in computer graphics, we assume light contains an equal amount of s-polarized and p-polarized light waves. Then the Fresnel reflection coefficient is simply an average between the s-polarized and p-polarized light reflection coefficients, as the comment indicates. If you have more information on how much s-polarized light is emitted by your light source, then the average should reflect that.

BRDF Explorer Sample

The following shader script implements the BRDF in the Disney BRDF Explorer tool, using the stock Blinn-Phong shader with the default microfacet distribution. Note how we just implemented the code separately and multiplied the BRDF by the modified "thin film" Fresnel term.

analytic

# Blinn Phong based on halfway-vector with single-layer thin
# film wave interference effects via a Fresnel film coating.

::begin parameters
float thickness 0 3000 250  # Thin film thickness (in nm)
float externalIOR 0.2 3 1     # External (air) refractive index
float thinfilmIOR 0.2 3 1.5   # Layer (thin film) refractive index
float internalIOR 0.2 3 1.25  # Internal (object) refractive index
float n 1 1000 100          # Blinn-Phong microfacet exponent
::end parameters

::begin shader

const float PI = 3.14159265f;

/* Amplitude reflection coefficient (s-polarized) */
float rs(float n1, float n2, float cosI, float cosT)
{
    return (n1 * cosI - n2 * cosT) / (n1 * cosI + n2 * cosT);
}

/* Amplitude reflection coefficient (p-polarized) */
float rp(float n1, float n2, float cosI, float cosT)
{
    return (n2 * cosI - n1 * cosT) / (n1 * cosT + n2 * cosI);
}

/* Amplitude transmission coefficient (s-polarized) */
float ts(float n1, float n2, float cosI, float cosT)
{
    return 2 * n1 * cosI / (n1 * cosI + n2 * cosT);
}

/* Amplitude transmission coefficient (p-polarized) */
float tp(float n1, float n2, float cosI, float cosT)
{
    return 2 * n1 * cosI / (n1 * cosT + n2 * cosI);
}

/* Pass the incident cosine. */
vec3 FresnelCoating(float cos0)
{
    /* Precompute the reflection phase changes (depends on IOR) */
    float delta10 = (thinfilmIOR < externalIOR) ? PI : 0.0f;
    float delta12 = (thinfilmIOR < internalIOR) ? PI : 0.0f;
    float delta = delta10 + delta12;

    /* Calculate the thin film layer (and transmitted) angle cosines. */
    float sin1 = pow(externalIOR / thinfilmIOR, 2) * (1 - pow(cos0, 2));
    float sin2 = pow(externalIOR / internalIOR, 2) * (1 - pow(cos0, 2));
    if ((sin1 > 1) || (sin2 > 1)) return vec3(1); /* Account for TIR. */
    float cos1 = sqrt(1 - sin1), cos2 = sqrt(1 - sin2);

    /* Calculate the interference phase change. */
    vec3 phi = vec3(2 * thinfilmIOR * thickness * cos1);
    phi *= 2 * PI / vec3(650, 510, 475);
    phi += delta;

    /* Obtain the various Fresnel amplitude coefficients. */
    float alpha_s = rs(thinfilmIOR, externalIOR, cos1, cos0)
                  * rs(thinfilmIOR, internalIOR, cos1, cos2);
    float alpha_p = rp(thinfilmIOR, externalIOR, cos1, cos0)
                  * rp(thinfilmIOR, internalIOR, cos1, cos2);
    float beta_s  = ts(externalIOR, thinfilmIOR, cos0, cos1)
                  * ts(thinfilmIOR, internalIOR, cos1, cos2);
    float beta_p  = tp(externalIOR, thinfilmIOR, cos0, cos1)
                  * tp(thinfilmIOR, internalIOR, cos1, cos2);

    /* Calculate the s- and p-polarized intensity transmission coefficient. */
    vec3 ts = pow(beta_s, 2) / (pow(alpha_s, 2) - 2 * alpha_s * cos(phi) + 1);
    vec3 tp = pow(beta_p, 2) / (pow(alpha_p, 2) - 2 * alpha_p * cos(phi) + 1);

    /* Calculate the transmitted power ratio for medium change. */
    float beamRatio = (internalIOR * cos2) / (externalIOR * cos0);

    /* Calculate the average reflectance. */
    return 1 - beamRatio * (ts + tp) * 0.5f;
}

vec3 BRDF(vec3 L, vec3 V, vec3 N, vec3 X, vec3 Y)
{
    vec3 H = normalize(L + V);
    float val = pow(max(0, dot(N, H)), n);
    return vec3(val) * FresnelCoating(dot(V, H));
}

::end shader

It is worth noting that this is a reference implementation meant to be readable, and can be thoroughly optimized. In particular, the Fresnel calculations are the most expensive, but there are numerous ways of reducing the amount of computations. For instance, we can use the reciprocity properties of s-polarized light, and also recycle many intermediate calculations. If you are not interested in perfect physical accuracy, you can also skip the polarization calculations and directly use intensity Fresnel coefficients, though because amplitudes are signed and intensities are not, you will need to calculate the proper sign to use for the cosine term somehow (or just ignore it altogether and have incorrect but plausible thin film interference).

If you are really desperate about runtime performance, you can still retain the nice colorful patterns while trading physical accuracy by approximating the final formula however you see fit, the only fundamental requirement is that the term be in there somewhere.

This is a screenshot of the above BRDF's polar plot at incidence 45 degrees and illustrates its wavelength-dependent nature:

Note how the BRDF differs for the three channels (in fact, every wavelength produces a different response, but we're working in RGB mode here). And here are a few renders (still from BRDF Explorer) with some sensible parameters. Here we assume the internal medium is fully opaque:

What happens when we set the thin film thickness to zero? In this case, the layer physically disappears and the formula degenerates to ordinary Fresnel reflection (more specifically, the Fresnel reflection coefficients for the layer become zero while the transmission coefficients become one).

What about making the thin film extremely thick? In that case, we see that the rate of change of the phase with respect to view angle becomes arbitrarily large, causing the interference effects to average out to white light, as expected. Also, because we are using a BRDF, we are assuming that light exits the surface at the same point it enters it, which is a very good approximation when the thin film is very small (on the order of light's wavelength). However, as the film becomes thicker, the approximation breaks down, so the film should probably be no larger than a few thousand nanometers, at most.

The same holds true for microfacet distributions. Thin films coated over surfaces with very high microfacet roughness coefficients are somewhat unphysical, since a coating naturally tends to be smoother than the surface it is applied on. This should be kept in mind, as the two layer interfaces are assumed to be coplanar.

You might also wonder what happens if refractive index depends on wavelength. Well, not much, the correct refractive index and incident/transmitted angles are simply used, and everything else remains the same, since waves of different frequencies do not interfere in any meaningful way. The BRDF above chooses to assume a constant IOR, though, to simplify matters. Also, if the refractive indices are wavelength-dependent, you will also observe dispersion effects in the transmitted light.

Transparency?

You may also want to handle transparency if the internal medium is not opaque. You can use whatever method you already have in place to render refractive surfaces, using the final transmitted angle (cos2 in the pseudocode). This is necessary for soap bubbles. Of course, if the internal medium is opaque, this is not necessary as the transmitted light is simply absorbed. It is also possible to use this with subsurface scattering (thus representing a subsurface scattering material with a thin film coating) by using the transmitted light (suitably refracted, as mentioned above).

Here is a render of a model with a soap-bubble-like BRDF, rendered with ray tracing. In this case, there is no visible refraction because soap bubbles are simply an air/water/air interface, so the final transmitted angle is the same as the incident angle:

Final notes

A good selection of parameters is essential to obtain realistic results. For instance, the film thickness should be on the order of light's wavelength (a few hundred nanometers). As you increase the thickness, interference effects disappear and as the thickness tends to zero, you just get ordinary Fresnel reflection, as mentioned previously. Make sure to use correct refractive indices for your materials. The range of values which can produce interference effects is quite narrow, so the parameters have to be accurate.

For metals or materials where the refractive index varies considerably over the visible spectrum, such as copper, three refractive indices (one per RGB channel) should be used for physical accuracy if possible. This requires only minor changes to the BRDF, as everything can be vectorized. It suffices to make the IOR parameters 3-component vectors and vectorize the Fresnel coefficient functions. The computational cost is therefore exactly the same.

One last point is that for non-solid thin films, such as oil or water coatings, the thickness of the layer is probably not constant at every point of a given object. As an example, soap bubbles are thicker at the bottom than at the top, due to gravity. For a convincing render, this should be taken into account. As a result, thin film thickness should probably be a vertex attribute rather than a material attribute, or, alternatively, a more general reflectance model should be considered (such as a spatially varying BRDF). Adding some noise to the film thickness can also go very far in improving the appearance of some materials, and it is convenient to implement.

Attached is the zipped BRDF Explorer script so that you may play around with it at your leisure.

Article Update Log

2 May 2013: Fixed a couple of bugs in the shader, corrected a few typos and improved formatting.
30 Apr 2013: Added some notes about interesting variations to apply to film thickness and on optimizations.
29 Apr 2013: Added notes on the motivation of neglecting absorption effects.
28 Apr 2013: Added notes on IOR and physical accuracy of solution.
27 Apr 2013: Added extra render, improved formatting.
26 Apr 2013: Added BRDF and some renders.
25 Apr 2013: Began writing article.

2D Animation Basics

Thu, 13 Dec 2012 05:59:29 +0000

Introduction

Simply put, 2D animation is movement and transformation of objects on the screen in two dimensions. A good example of 2D animation is classic cartoons – multiple pictures of Mickey Mouse or Donald Duck alternating over time and producing the effect of moving objects. The modern 2D animation software allows creating animation much easier though, without making oodles of frames.

Frames and Keyframes

Just like with classic animation, computer 2D animation also consists of frames. Frame is a single picture of an object illustrating certain position, size and other properties of that object. Adjacent frames are slightly different from each other, so when the animation is played at the full speed – usually 24 frames per second – the effect of motion appears. This is similar to cartoons and movies.

But in computer animation you typically don’t have to draw each and every frame manually. Most of 2D animation software allows you to set up keyframes – frames that specify intermediate positions of an object. The motion of the object between these keyframes is built by the software automatically. So you end up with creating keyframes only, while the animation tool calculates and builds the rest of frames automatically. This process is called interpolation.

Curve Animation

The motion between two keyframes can be different – accelerating, slowing down, steady or even all of these in one move. The simplest way to specify how exactly one keyframe should translate to another is curves. A curve defines how a specific parameter of an object, such as a coordinate or a scale should vary over time.

Let’s take a look at a simple animation example made in one of 2D animation programs – GameDev Animation Studio. We want to make a ball move from up to down and then bounce back.

Basically, this means the vertical coordinate of the object, or simply the Y-coordinate should change from zero (the top of the screen) to the maximum value that corresponds to the bottom of the screen. Then the ball “bounces”, i.e. its form gets distorted in the vertical direction and restored again. After that, the ball moves up, which means the Y-coordinate should decrease back to zero. Also note, that due to gravity, the ball moves down with acceleration, while the backward movement must be decelerating. This is mirrored in the slope of the curve: a sharper slope means faster movement and vise versa.

Position Curve

So the curve describing vertical movement of the ball looks as follows:

As you see, the curve here gradually develops from zero to its maximum, and then, after the peak, it goes back to zero again. The highest value of the curve corresponds to the moment the ball hits the ground – this happens in the 30th frame on the timeline. After that moment, the ball starts moving up.

Zoom Curve

When the ball is just to hit the ground, it starts deforming. Indeed, the ball is elastic, so it kinda must deform. This animation can be done with another curve – the zoom curve. We’ll change the vertical scale (Y-scale) of the ball to reflect the way it gets distorted during the hit.

Above is the Y-scale zoom curve. As you see, prior the 29th frame the curve is at its maximum. This corresponds to the non-deformed ball falling down. In 29th, the curve goes down, and the ball’s image is zoomed down along the vertical axis. Then, after the curve passes the “hit-the-ground” point in 31st frame, the image is zoomed up to normal, so the curve goes up again to its maximum.

Note, that the higher the zoom factor is, the more elastic the ball will look. Thus, experimenting with curves in 2D animation tools is an easy way to find the most spectacular and satisfying animation effects. On the other hand, the classic animation approach would require you to discard all the work already done and start from scratch if you decided to change something.

Adding Details to Animation

Obviously, the ball just moving up and down doesn’t look too realistic, that is why we need a shadow. We should place a shadow picture beneath the ball and use curves to make it look real. Overall, the process is pretty much similar to above steps. We’ll need to scale the shadow and modify its transparency level in accordance with the ball’s movement.

Here is the final result of the jumping ball animation.

Conclusion

Making 2D animation is a piece of cake with proper software. Thanks to automated image transformation between keyframes, curve animation and multiple transition effects, creating an animated picture or a movie can be done in several minutes, even by a beginner.

GPGPU image processing basics using OpenCL.NET

Thu, 29 Nov 2012 10:09:42 +0000

OpenCL is a cross-platform framework used mostly for GPGPU (General-purpose computing on graphics processing units). There are plenty of tutorials available on image processing with OpenCL using C/C++, however there's not much information that would cover OpenCL image processing with .NET.

I won't go into details about OpenCL kernels/queues/etc. (there's plenty of information available on the internet), however I'll provide you with a bare minimum code required to load an image from disk, process it with OpenCL on the GPU and save it back to a file.

Before we get started, make sure that you download the source code of OpenCL.NET from http://openclnet.codeplex.com/ and add it to your project.

We'll use a simple OpenCL kernel that converts an input image into a grayscale image. The kernel should be saved to a separate file.

Kernel source code:

__kernel void imagingTest(__read_only  image2d_t srcImg,
					   __write_only image2d_t dstImg)
{
  const sampler_t smp = CLK_NORMALIZED_COORDS_FALSE | //Natural coordinates
	CLK_ADDRESS_CLAMP_TO_EDGE | //Clamp to zeros
	CLK_FILTER_LINEAR;
  int2 coord = (int2)(get_global_id(0), get_global_id(1));

  uint4 bgra = read_imageui(srcImg, smp, coord); //The byte order is BGRA

  float4 bgrafloat = convert_float4(bgra) / 255.0f; //Convert to normalized [0..1] float

  //Convert RGB to luminance (make the image grayscale).
  float luminance =  sqrt(0.241f * bgrafloat.z * bgrafloat.z + 0.691f * bgrafloat.y * bgrafloat.y + 0.068f * bgrafloat.x * bgrafloat.x);
  bgra.x = bgra.y = bgra.z = (uint) (luminance * 255.0f);

  bgra.w = 255;

  write_imageui(dstImg, coord, bgra);
}

Namespaces Used

using System;
using System.Collections;
using System.Collections.Generic;
using System.Drawing;
using System.Drawing.Imaging;
using System.IO;
using System.Runtime.InteropServices;
using OpenCL.Net;

Error handling

Since OpenCL.NET is a wrapper for C API, we'll have to do all the error checking on our own. I'm using the following two methods:

private void CheckErr(Cl.ErrorCode err, string name)
{
	if (err != Cl.ErrorCode.Success) {
		Console.WriteLine("ERROR: " + name + " (" + err.ToString() + ")");
	}
}

private void ContextNotify(string errInfo, byte[] data, IntPtr cb, IntPtr userData) {
	Console.WriteLine("OpenCL Notification: " + errInfo);
}

Setting Up

The following two variables should be declared in the class itself and will be shared across all of the methods:

private Cl.Context _context;
private Cl.Device _device;

And this is the method that sets up OpenCL:

private void Setup ()
{
	Cl.ErrorCode error;
	Cl.Platform[] platforms = Cl.GetPlatformIDs (out error);
	List devicesList = new List ();
  
	CheckErr (error, "Cl.GetPlatformIDs");
  
	foreach (Cl.Platform platform in platforms) {
		string platformName = Cl.GetPlatformInfo (platform, Cl.PlatformInfo.Name, out error).ToString ();
		Console.WriteLine ("Platform: " + platformName);
		CheckErr (error, "Cl.GetPlatformInfo");

		//We will be looking only for GPU devices
		foreach (Cl.Device device in Cl.GetDeviceIDs(platform, Cl.DeviceType.Gpu, out error)) {
			CheckErr (error, "Cl.GetDeviceIDs");
			Console.WriteLine ("Device: " + device.ToString ());
			devicesList.Add (device);
		}
	}
  
	if (devicesList.Count <= 0) {
		Console.WriteLine ("No devices found.");
		return;
	}
  
	_device = devicesList[0];
  
	if (Cl.GetDeviceInfo(_device, Cl.DeviceInfo.ImageSupport, out error).CastTo() == Cl.Bool.False)
	{
		Console.WriteLine("No image support.");
		return;
	}

	_context = Cl.CreateContext(null, 1, new[] { _device }, ContextNotify, IntPtr.Zero, out error);	//Second parameter is amount of devices
	CheckErr(error, "Cl.CreateContext");
}

The Image Processing Part

The main problem is that OpenCL.NET is a wrapper around C API of OpenCL, so it can only work with unmanaged memory. However, all of the data in .NET is managed, so we'll have to marshal the data between managed/unmanaged memory. Usually it would be much easier to handle the RGBA color components in float [0..1] space. However, the input image should be in byte[] array, because it would really affect the performance to do the byte=>float conversion on the CPU (we would have to divide each component by 255 for every pixel of the image twice - once before the image processing and once after).

public void ImagingTest (string inputImagePath, string outputImagePath)
{
	Cl.ErrorCode error;

	//Load and compile kernel source code.
	string programPath = Environment.CurrentDirectory + "/../../ImagingTest.cl";	//The path to the source file may vary
  
	if (!System.IO.File.Exists (programPath)) {
		Console.WriteLine ("Program doesn't exist at path " + programPath);
		return;
	}
  
	string programSource = System.IO.File.ReadAllText (programPath);
  
	using (Cl.Program program = Cl.CreateProgramWithSource(_context, 1, new[] { programSource }, null, out error)) {
		CheckErr(error, "Cl.CreateProgramWithSource");

		//Compile kernel source
		error = Cl.BuildProgram (program, 1, new[] { _device }, string.Empty, null, IntPtr.Zero);
		CheckErr(error, "Cl.BuildProgram");

		//Check for any compilation errors
		if (Cl.GetProgramBuildInfo (program, _device, Cl.ProgramBuildInfo.Status, out error).CastTo()
			!= Cl.BuildStatus.Success) {
			CheckErr(error, "Cl.GetProgramBuildInfo");
			Console.WriteLine("Cl.GetProgramBuildInfo != Success");
			Console.WriteLine(Cl.GetProgramBuildInfo(program, _device, Cl.ProgramBuildInfo.Log, out error));
			return;
		}

		//Create the required kernel (entry function)
		Cl.Kernel kernel = Cl.CreateKernel(program, "imagingTest", out error);
		CheckErr(error, "Cl.CreateKernel");
	  
		int intPtrSize = 0;
		intPtrSize = Marshal.SizeOf(typeof(IntPtr));

		//Image's RGBA data converted to an unmanaged[] array
		byte[] inputByteArray;
		//OpenCL memory buffer that will keep our image's byte[] data.
		Cl.Mem inputImage2DBuffer;

		Cl.ImageFormat clImageFormat = new Cl.ImageFormat(Cl.ChannelOrder.RGBA, Cl.ChannelType.Unsigned_Int8);

		int inputImgWidth, inputImgHeight;
	  
		int inputImgBytesSize;

		int inputImgStride;

		//Try loading the input image
		using (FileStream imageFileStream = new FileStream(inputImagePath, FileMode.Open) ) {
			System.Drawing.Image inputImage = System.Drawing.Image.FromStream( imageFileStream );
		  
			if (inputImage == null) {
				Console.WriteLine("Unable to load input image");
				return;
			}
		  
			inputImgWidth = inputImage.Width;
			inputImgHeight = inputImage.Height;
		  
			System.Drawing.Bitmap bmpImage = new System.Drawing.Bitmap(inputImage);

			//Get raw pixel data of the bitmap
			//The format should match the format of clImageFormat
			BitmapData bitmapData = bmpImage.LockBits( new Rectangle(0, 0, bmpImage.Width, bmpImage.Height),
													  ImageLockMode.ReadOnly, PixelFormat.Format32bppArgb);//inputImage.PixelFormat);

			inputImgStride = bitmapData.Stride;
			inputImgBytesSize = bitmapData.Stride * bitmapData.Height;
		  
			//Copy the raw bitmap data to an unmanaged byte[] array
			inputByteArray = new byte[inputImgBytesSize];
			Marshal.Copy(bitmapData.Scan0, inputByteArray, 0, inputImgBytesSize);

			//Allocate OpenCL image memory buffer
			inputImage2DBuffer = Cl.CreateImage2D(_context, Cl.MemFlags.CopyHostPtr | Cl.MemFlags.ReadOnly, clImageFormat,
												(IntPtr)bitmapData.Width, (IntPtr)bitmapData.Height,
												(IntPtr)0, inputByteArray, out error);
			CheckErr(error, "Cl.CreateImage2D input");
		}

		//Unmanaged output image's raw RGBA byte[] array
		byte[] outputByteArray = new byte[inputImgBytesSize];

		//Allocate OpenCL image memory buffer
		Cl.Mem outputImage2DBuffer = Cl.CreateImage2D(_context, Cl.MemFlags.CopyHostPtr | Cl.MemFlags.WriteOnly, clImageFormat,
													  (IntPtr)inputImgWidth, (IntPtr)inputImgHeight, (IntPtr)0, outputByteArray, out error);
		CheckErr(error, "Cl.CreateImage2D output");

		//Pass the memory buffers to our kernel function
		error = Cl.SetKernelArg(kernel, 0, (IntPtr)intPtrSize, inputImage2DBuffer);
		error |= Cl.SetKernelArg(kernel, 1, (IntPtr)intPtrSize, outputImage2DBuffer);
		CheckErr(error, "Cl.SetKernelArg");
	  
		//Create a command queue, where all of the commands for execution will be added
		Cl.CommandQueue cmdQueue = Cl.CreateCommandQueue(_context, _device, (Cl.CommandQueueProperties)0, out error);
		CheckErr(error, "Cl.CreateCommandQueue");

		Cl.Event clevent;

		//Copy input image from the host to the GPU.
		IntPtr[] originPtr = new IntPtr[] { (IntPtr)0, (IntPtr)0, (IntPtr)0 };	//x, y, z
		IntPtr[] regionPtr = new IntPtr[] { (IntPtr)inputImgWidth, (IntPtr)inputImgHeight, (IntPtr)1 };	//x, y, z
		IntPtr[] workGroupSizePtr = new IntPtr[] { (IntPtr)inputImgWidth, (IntPtr)inputImgHeight, (IntPtr)1 };
		error = Cl.EnqueueWriteImage(cmdQueue, inputImage2DBuffer, Cl.Bool.True, originPtr, regionPtr, (IntPtr)0, (IntPtr)0, inputByteArray, 0, null, out clevent);
		CheckErr(error, "Cl.EnqueueWriteImage");

		//Execute our kernel (OpenCL code)
		error = Cl.EnqueueNDRangeKernel(cmdQueue, kernel, 2, null, workGroupSizePtr, null, 0, null, out clevent);
		CheckErr(error, "Cl.EnqueueNDRangeKernel");

		//Wait for completion of all calculations on the GPU.
		error = Cl.Finish(cmdQueue);
		CheckErr(error, "Cl.Finish");

		//Read the processed image from GPU to raw RGBA data byte[] array
		error = Cl.EnqueueReadImage(cmdQueue, outputImage2DBuffer, Cl.Bool.True, originPtr, regionPtr,
									(IntPtr)0, (IntPtr)0, outputByteArray, 0, null, out clevent);
		CheckErr(error, "Cl.clEnqueueReadImage");

		//Clean up memory
		Cl.ReleaseKernel(kernel);
		Cl.ReleaseCommandQueue(cmdQueue);
	  
		Cl.ReleaseMemObject(inputImage2DBuffer);
		Cl.ReleaseMemObject(outputImage2DBuffer);

		//Get a pointer to our unmanaged output byte[] array
		GCHandle pinnedOutputArray = GCHandle.Alloc(outputByteArray, GCHandleType.Pinned);
		IntPtr outputBmpPointer = pinnedOutputArray.AddrOfPinnedObject();

		//Create a new bitmap with processed data and save it to a file.
		Bitmap outputBitmap = new Bitmap(inputImgWidth, inputImgHeight, inputImgStride, PixelFormat.Format32bppArgb, outputBmpPointer);
	  
		outputBitmap.Save(outputImagePath, System.Drawing.Imaging.ImageFormat.Png);

		pinnedOutputArray.Free();
	}
}

Now you should have a good foundation for more complex image processing effects on the GPU.

Dynamic Resolution Rendering

Mon, 26 Sep 2011 19:35:42 +0000

Brought to you by the Intel® Visual Computing Developer Community | Download the source code | Watch the video

Introduction

The resolution selection screen has been one of the defining aspects of PC gaming since the birth of 3D games. In this whitepaper and the accompanying sample code, we argue that this no longer needs to be the case; developers can dynamically vary the resolution of their rendering instead of having a static resolution selection.

Dynamic resolution rendering involves adjusting the resolution to which you render the 3D scene by constraining the rendering to a portion of a render target using a viewport, and then scaling this to the output back buffer. Graphical user interface components can then be rendered at the back buffer resolution, as these are typically less expensive elements to draw. The end result is that stable high frame rates can be achieved with high quality GUIs.

We'll be presenting performance results and screenshots in this article taken on a pre-release mobile 2nd generation Intel® Core™ i7 processor (Intel® microarchitecture code name Sandy Bridge, D1 stepping quad core 2.4 GHz CPU with 4GB DDR3 1333MHz RAM) with Intel® HD Graphics 3000.

This article and the accompanying sample were originally presented at the Game Developers Conference (GDC) in San Francisco 2011, and a video of the presentation can be found on GDC Vault [GDC Vault 2011], with the slides for that presentation available on the Intel website [Intel GDC 2011]. Since the presentation, the author has discovered that several game companies already use this technique on consoles; Dmitry Andreev from LucasArts' presentation on Anti-Aliasing is the only public source, though with few details on the dynamic resolution technique used [Andreev 2011].

Figure 1: The sample scene viewed from one of the static camera viewpoints.

  Motivation

  Games have almost always had a strong performance variation with resolution, and the increase in shader complexity along with post-processing techniques has continued the trend of per-pixel costs dominating modern games. Increasing resolution also increases texture sampling and render target bandwidth. Setting the resolution appropriately for the performance of the system is therefore critical. Being able to vary the resolution dynamically gives the developer an additional performance control option which can enable the game to maintain a stable and appropriate frame rate, thus improving the overall quality of the experience.

Rendering the graphical user interface at the native screen resolution can be particularly important for role playing, real time strategy, and massively multiplayer games. Suddenly, even on low-end systems, the player can indulge in complex chat messaging whilst keeping an eye on their teammates' stats.

Finally, with the increasing dominance of laptops in PC gaming, power consumption is beginning to become relevant to game development. Performance settings can cause a reduction in CPU and GPU frequency when a machine goes from mains to battery power, and with dynamic resolution rendering, the game can automatically adjust the resolution to compensate. Some games may want to give the user the option of a low power profile to further reduce power consumption and enable longer gaming on the go. Experiments with the sample have found that cutting the resolution to 0.5x reduces the power consumption of the processor package to 0.7x normal when vertical sync is enabled so that the frame rate is maintained.

Basic Principles

  The basic principle of dynamic resolution rendering is to use a viewport to constrain the rendering to a portion of an off-screen render target, and then to scale the view. For example, the render target might be of size (1920, 1080), but the viewport could have an origin of (0, 0) and size (1280, 720).

Figure 2: using a viewport to constrain rendering

By creating render targets larger than the back buffer, the dynamic resolution can be varied from subsampled to supersampled. Care needs to be taken to ensure the full set of required render targets and textures fit within graphics memory, but systems based on Intel® microarchitecture code name Sandy Bridge processor graphics usually have considerable memory, as they use system memory.

Figure 3: dynamic resolution can be varied from subsampling to supersampling

When undertaking normal rendering to the dynamic viewport, there are no changes that need to be made-the rasterization rules ensure this is handled. However, when reading from the render target, care needs to be taken to scale the coordinates appropriately and handle clamping at the right and bottom edges.

The following example pixel shader code shows how to clamp UVs. This is mainly used when doing dependent reads (i.e., when there are per-pixel operations on a UV, which is subsequently used to sample from a dynamic render target).

In the case of motion blur-a common post-process operation that uses dependent reads from a render target-the extra math required has little effect on the performance, as the shader is texture-fetch bound.

Figure 4: Color leak on edges of screen due to motion blur, which can be solved by using clamping

In addition to clamping, it's also important to ensure that the resolution ratios used in shaders is representative of the actual viewport ratio, rather than just your application's desired ratio. This is easily obtained by recalculating the ratio from the dynamic viewport dimensions. For example, in the sample code function DynamicResolution::SetScale, the following is performed after ensuring the scale meets boundary criteria:

Scaling Filters
After rendering the 3D scene, the viewport area needs to be scaled to the back buffer resolution. A variety of filters can be used to perform this, and the sample implements several examples as described here.

Point Filtering
Point filtering is a fast basic filter option. Scaling from a 0.71x ratio dynamic viewport to 1280x720 takes ~0.4ms.

Bilinear Filtering
Bilinear filtering is almost as fast as point filtering due to hardware support, and it reduces the aliasing artifacts from edges by smoothing, but also blurs the scene. Scaling from a 0.71x ratio dynamic viewport to 1280x720 takes ~0.4ms.

Bicubic Filtering
Bicubic filtering is only noticeably better than bilinear for resolutions of 0.5x the back buffer, and its performance is 7x slower even using a fast bicubic filter [Sigg 2005]. Scaling from a 0.71x ratio dynamic viewport to 1280x720 takes ~3.5ms.

Noise Filtering
Adding some noise to point filtering helps to add high frequencies, which break the aliasing slightly at a low cost. The implementation in the sample is fairly basic, and improved film grain filtering might artistically fit your rendering. Scaling from a 0.71x ratio dynamic viewport to 1280x720 takes ~0.5ms.

Noise Offset Filtering
Adding a small random offset to the sampling location during scaling reduces the regularity of aliased edges. This approach is common in fast filtering of shadow maps. Scaling from a 0.71x ratio dynamic viewport to 1280x720 takes ~0.7ms.

Temporal Anti-aliasing Filtering
This scaling filter requires extra support during the initial rendering path to render odd and even frames offset by half a pixel in X and Y. When filtered intelligently to remove ghosting artifacts, the resulting image quality is substantially improved by sampling from twice as many pixels. This filtering method is described in greater depth in its own section below. Scaling from a 0.71x ratio dynamic viewport to 1280x720 takes ~1.1ms, and has almost the same quality as rendering to full resolution.

Temporal Anti-aliasing Details
Temporal anti-aliasing has been around for some time; however, ghosting problems due to differences in the positions of objects in consecutive frames have limited its use. Modern rendering techniques are finally making it an attractive option due to its low performance overhead.

The basic approach is to render odd and even frames jittered (offset) by half a pixel in both X and Y. The sample code does this by translating the projection matrix. The final scaling then combines both the current and previous frames, offsetting them by the inverse of the amount they were jittered. The final image is thus made from twice the number of pixels arranged in a pattern similar to the dots of the five side on a die, frequently termed a quincunx pattern.

Figure 5: Temporal Anti-Aliasing basic principle

Used along with dynamic resolution, this approach gives an increased observed number of pixels in the scene when the dynamic resolution is lower than the back buffer, improving the detail in the scene. When the dynamic resolution is equal or higher to the back buffer, the result is a form of anti-aliasing.

Figure 6: Result of Temporal AA when dynamic resolution is lower than that of the back buffer

Figure 7: Result of Temporal AA when dynamic resolution is equal or higher to that of the back buffer

In order to get increased texture resolution, a MIP LOD bias needs to be applied to textures. In Microsoft Direct3D* 11, use a D3D11_SAMPLER_DESC MipLODBias of -0.5f during the 3D scene pass. Additionally, the sampler used during scaling needs to use bilinear minification filtering, for example: D3D11_FILTER_MIN_LINEAR_MAG_MIP_POINT.

In order to reduce ghosting, we use the velocity buffer written out for motion blur. Importantly, this buffer contains the velocity for each pixel in screen space, thus accounting for camera movement. A scale factor is calculated from both the current and previous frame's velocity and applied to the previous frame's colour to determine its contribution to the final image. This scales the contribution based on how similar the sample location is in real space in both frames.

The sample has K tuned to give what the author considers to be the best results for a real time application, with no ghosting observed at realistically playable frame rates. Screenshots do expose a small amount of ghosting in high contrast areas as in the screenshot below, which can be tuned out if desired.

For games, transparencies present a particular problem in not always rendering out velocity information. In this case, the alpha channel could be used during the forwards rendering of the transparencies to store a value used to scale the contributions in much the same way as the velocity is currently used.

An alternative to this approach for ghosting removal is to use the screen space velocity to sample from the previous frame at the location where the current pixel was. This is the technique used in the CryENGINE* 3, first demonstrated in the game Crysis* 2 [Crytek 2010]. Intriguingly, LucasArts' Dmitry Andreev considered using temporal anti-aliasing, but did not due to the use of dynamic resolution in their engine [Andreev 2011]. The author believes these are compatible, as demonstrated in the sample code.

Figure 8: Temporal Anti-Aliasing with velocity scaling and moving objects

  The Effect of Motion Blur

Motion blur smears pixels and reduces observed aliasing effectively, hence a lower resolution can be used when the camera is moving. However, the sample does not exploit this in its resolution control algorithm. The following screenshots show how reducing the resolution to 0.71x the back buffer results in higher performance, but roughly the same image. Combined with varying motion blur sample rates, this could be a way to reduce artifacts from undersampling with large camera motions whilst maintaining a consistent performance.

Figure 9: Motion blur with dynamic resolution off

Figure 10: Motion blur with dynamic resolution on at 0.71x resolution. Note the decreased frame time yet similar quality end result

  Supersampling

  Supersampling is a simple technique where the render target used to render the scene is larger than the back buffer. This technique is largely ignored by the current real-time rendering community-multi sampled anti-aliasing and other anti-aliasing techniques have replaced its usage due to their better memory consumption and performance.

Using dynamic resolution significantly reduces the performance impact of adding supersampling, as the actual resolution used can be dynamically adjusted. There is a small performance impact to enabling supersampling, mainly due to the extra cost of clearing the larger buffers. The sample code implements a 2x resolution render target when supersampling is enabled, but good quality results are observed for relatively small increases in resolution over the back buffer resolution, so a smaller render target could be used if memory were at a premium. Memory is less of an issue on processor graphics platforms, as the GPU has access to a relatively large proportion of the system memory, all of which is accessible at full performance.

Once dynamic resolution rendering methods are integrated, using supersampling is trivial. We encourage developers to consider this, since it can be beneficial for smaller screen sizes and future hardware which could have sufficient performance to run the game at more than its maximum quality.

  Render Target Clearing

  Since dynamic resolution rendering does not always use the entire render targets surface, it can be beneficial to clear only the required portion. The sample implements a pixel shader clear, and on the Intel® HD Graphics 3000-based system tested, the performance of a pixel shader clear was greater than that of a standard clear when the dynamic ratio was less than 0.71x for a 1280x720 back buffer. In many cases, it may not be necessary to clear the render targets, as these get overwritten fully every frame.

Depth buffers should still be cleared completely with the standard clear methods, since these may implement hierarchical depth. Some multi-sampled render targets may also use compression, so should be cleared normally.

  Performance Scaling

  The sample code scales well with resolution, despite the heavy vertex processing load due to the large highly detailed scene with no level of detail and only very simple culling performed. This gives the chosen control method significant leverage to maintain frame rate at the desired level.

Most games use level-of-detail mechanisms to control the vertex load. If these are linked to the approximate size of the object in pixels, the resulting performance scaling will be greater.

Figure 11: Dynamic Resolution Performance at 1280x720

  Resolution Control

  The sample implements a resolution control method in addition to allowing manual control. The code is in the file DynamicResolutionRendering.cpp, in the function ControlResolution. The desired performance can be selected between the refresh rate (usually 60Hz or 60FPS) and half the refresh rate (usually 30FPS).

The control scheme is basic: a resolution scale delta is calculated proportionally to the dimensionless difference in the desired frame time and the current frame time.

Where S' is the new resolution scale ratio, S is the current resolution scale ratio, is the scale delta, k a rate of change constant, T the desired frame time, and t the current frame time.

The current frame time uses an average of the GPU inner frame time excluding the present calculated using Microsoft DirectX* queries, and the frame time calculated from the interval between frames in the normal way. The GPU inner frame time is required when vertical sync is enabled, as in this situation the frame time is capped to the sync rate, yet we need to know if the actual rendering time is shorter than that. Averaging with the actual frame rate helps to take into account the present along with some CPU frame workloads. If the actual frame time is significantly larger than the GPU inner frame time, this is ignored, as these are usually due to CPU side spikes such as going from windowed to fullscreen.

Potential Improvements

The following list is by no means complete, but merely some of the features which the author believes would naturally extend the current work:

Combine the dynamic resolution scene rendering with a similar method for shadow maps.
Use this technique with a separate control mechanism for particle systems, allowing enhanced quality when only a few small particles are being rendered and improved performance when the fill rate increases.
The technique is compatible with other anti-aliasing techniques that can also be applied along with temporal anti-aliasing.
Temporal anti-aliasing can use an improved weighted sum dependent on the distance to the pixel center of the current and previous frames, rather than just a summed blend. A velocity-dependent offset read, such as that used in the CryENGINE* 3 [Crytek 2010], could also be used.
Some games may benefit from running higher quality anti-aliasing techniques over a smaller area of the image, such as for the main character or on RTS units highlighted by the mouse.

Conclusion

Dynamic resolution rendering gives developers the tools needed to improve overall quality with minimal user intervention, especially when combined with temporal anti-aliasing. Given the large range of performance in the PC GPU market, we encourage developers to use this technique as one of their methods of achieving the desired frame rate for their game.

References

[Sigg 2005] Christian Sigg, Martin Hadwiger, "Fast Third Order Filtering", GPU Gems 2. Addison-Wesley, 2005.

[Crytek 2010] HPG 2010 "Future graphics in games", Cevat Yerli & Anton Kaplanyan. http://www.crytek.co...e/presentations

[GDC Vault 2011] http://www.gdcvault....ution-Rendering

[Intel GDC 2011] http://software.inte...s/intelgdc2011/

[Andreev 2011] http://www.gdcvault....rom-a-Different [PPT 4.6MB]

Accelerating iray with the NVIDIA Quadro 5000

Mon, 15 Aug 2011 19:37:57 +0000

Introduction

First of all, I'm a software guy. I use software, I write about software and I sometimes even dream about software, but I'm also not stupid. So, when NVIDIA contacts me and says we've got a graphics board that makes your software a lot faster and we'd like you to try it out, my answer is a quick ‘yes.’

  The board they are referring to is the NVIDIA Quadro 5000. This is a professional level board with 2.5 GB of memory, and enabled with CUDA. It is this technology that is responsible for speeding up software including 3ds Max, Maya and many of the products in Adobe's Creative Suite, especially the Mercury Playback Engine contained in Premiere Pro.

  The Installation

  As a digital designer, it is my job to make images and animations look great using the amazing software at my disposal. This is made possible using the large black box that sits faithfully at my side. This black box doesn't get my attention at all unless it occasionally starts whirring when a fan turns on or as a place to support my important papers. Needless to say, I'd be more interested in alphabetizing my CD collection than cracking open my computer to install a new graphics card. But the promise of not only faster renders, but much faster renders, makes it worth the hassle.

  The NVIDIA Quadro 5000 graphics card, shown in Figure 1, is not for the faint of heart. It is a sizable card taking up two PCI slots and the entire card is encased like a prototype car coming out of Detroit that screams, "serious stuff is found under here." I was able to get the card to fit, albeit tightly, in my standard Dell Studio XPS box, but I had to remove the hard drive in order to get the monster placed.

Figure 1: The NVIDIA Quadro 5000 graphics board is a professional level piece of hardware and overkill for gamers.

  Once the graphics board was installed, my machine booted right up without any trouble and after a quick download of the software drivers from NVIDIA's web site, I found my machine was working just as before the installation ordeal. With business as usual set in, I started about my work keeping in the back of my mind the fact that I now had a monster GPU ready to show its stuff.


  First Impressions

  Having the killer graphics board installed and working without messing up my system relieved my initial hardware fears, so my next thought was, "let's see what this thing can do." I immediately opened the latest version of 3ds Max and loaded a recent NVIDIA iray scene that I created.

  Before installing the Quadro 5000 card, I rendered this sample scene by simply setting iray to do an Unlimited render before retiring for the evening. I remember checking on the scene after an hour or so and it looked quite good, but I could still see a few artifacts.

  I rendered the same sample scene, shown in Figure 2, after installing the Quadro 5000 and I could tell it was cranking on the scene because I noticed a delightfully pleasant humming coming from my computer internals. I also noticed that the same quality of the final render was achieved while I watched in a matter of minutes instead of hours like before using the Quadro card. While I realize this isn't the typical quantitative benchmark that gives precise percentages, such as 10X faster; for me it means that I can render a scene without having to leave the computer.

Figure 2: This sample NVIDIA iray scene rendered to final quality in minutes rather than hours using the NVIDIA Quadro 5000.

  Tackling a Monster Project

  After testing out a relatively simple project, I decided to throw something with a little more weight at the new board. For my second render, I loaded a kitchen scene with lots of shiny, reflective surfaces made of stainless steel, porcelain and glass with lots of lights bouncing all over the place. The scene is actually an animation rendered in high-dev. For this test, I decided to set iray to do 100 passes per frame to see what kind of quality I'd get in a reasonable render time.

  The most complex frames of the animation took just over 7 minutes to complete the full 100 frames and the quality was a fair approximation of the ambient occlusion and all reflective surfaces, but the overall image was still grainy. After re-rendering the same frame set to 500 passes, the image was fairly close to photo-realistic with only a little graininess in the deep corners, and it took 44 minutes to render. Still a huge improvement from the earlier project render time, which took a over 24 hours for a 15 second animation.

  After several more test renders of my personal work, I decided to try out some of the scenes that NVIDIA provided me with, including a detailed Bugatti Veyron sports car, shown in Figure 3. The iray renderer for this car was set to do 500 passes at a resolution of 1300 by 900. The resulting image, shown in Figure 3, took just over 5 minutes. It also provided an excuse for me to include such a gorgeous image in this article.

Figure 3: The photo-realistic image quality of NVIDIA iray with the Quadro 5000 graphics card is amazing for taking only 5 minutes of render time.

  Understanding NVIDIA iray

  The speed and power behind using the Quadro board to render iray scenes is based on the fact that the Quadro board has been optimized specifically to work with iray by offloading much of the rendering process to the board's GPU. But, first we need to understand what iray is.

  When 3ds Max was created, it was purposefully made from an open architecture that allowed specific modules to be replaced by other 3rd party plug-ins. This allowed a young company, known as mental images, to do their magic by replacing the default rendering engine used by 3ds Max with their own rendering engine called mental ray (spelled without capital letters, just like the company, which was acquired by NVIDIA in 2008). Mental ray used many advanced rendering techniques to improve the render quality dramatically, especially for realistic scenes with lots of reflections. Over time, mental ray has pushed the limits of rendering quality and speed and is available for many different rendering packages and platforms.

  The one problem with mental ray is that with all its advanced settings, it can take a long time to figure out which settings give the best quality in the shortest amount of time. A problem that was exacerbated as the product became more and more complex. Mental images solution was to develop iray.

  NVIDIA iray is frightfully simple. It doesn't include all the advanced settings found in mental ray, but lets you set the render time or the number of passes and it automatically produces the best possible image given those constraints. It also has an unlimited option that will continually refine the render until it is stopped by the user. The one downside to iray is that it will only work with a unique set of materials that are available in 3ds Max. If a default material is used, it is simply rendered as flat white.

  3ds Max is a product developed by a team at Autodesk, and mental ray and iray is developed independent of Autodesk by NVIDIA. As mentioned, mental images was acquired in 2008 by a much bigger fish named NVIDIA, who has since put their team in contact with the brains behind the Quadro graphics card. The result is that it has been developed to take advantage of the latest features built into the Quadro graphics cards, enabling the software to render scenes with amazing quality faster than ever before and faster than any other graphics card.

Figure 4: This robot arm includes over 210,000 polys and was rendered in NVIDIA iray in just over 7 minutes.

  The Magic Behind the Technology

  So how is it done? The Quadro 5000 has a unit called the GigaThread Engine that lets several tasks be worked on at the same time including loading data into and out of the GPU while simultaneously doing some rendering processes. There is also a Scalable Geometry Engine that lets each cluster create triangles. By having each cluster work independently, the total throughput gets multiplied resulting in faster renders.

  The Quadro series of professional graphic cards also includes a feature called Error Correcting Code (ECC) that detects many multi-bit errors and re-runs the data to eliminate the errors. This is especially important in some critical areas like engineering and medical imaging. This feature can be disabled if the error correction isn't needed to make the process even faster.

  The Quadro 5000 has 352 Compute Unified Device Architecture (CUDA) parallel processor cores. It is the NVIDIA CUDA architecture that iray takes advantage of to speed its rendering cycles. The Quadro 5000 board also supports OpenGL 4.0, Shader Model 5.0 and DirectX 11. It is capable or rendering 950 million triangles per second.

  Other Benefits and Features

  The former mental images team, whose members are now are part of NVIDIA, isn't the only software product development team that has been working closely with NVIDIA engineers. The development team behind the Mercury Playback Engine found in the latest release of Adobe's Premiere Pro has also taken advantage of the optimized CUDA features. Quadro graphics boards allows video professionals to create and edit a large number of effects and transitions including color correction, blur, deinterlacing, compositing and blending, and play them back in real-time, which is a huge time-saver for video productions.

  The NVIDIA GPU found in the Quadro board also can accelerate physics simulations via the NVIDIA PhysX engine. 3ds Max's MassFX system includes an option to enable Hardware Acceleration that offloads some of the physics processing to the GPU.

  The Quadro 5000 graphics board also includes support for stereo 3D using the OpenGL stereo API. The NVIDIA Control Panel lets you set up and test out several different types of 3D glasses including the standard red/blue (anaglyph) glasses, and top-of-the-line active shutter 3D glasses like NVIDIA 3D Vision and 3D Vision Pro, the latter two of which NVIDIA recommends for viewing 3D games and applications. There is also a list of compatible games that you can play in 3D and a slide show of 3D game images.


  Summary

  In summary, I'm not at all surprised that the Quadro 5000 graphics board was faster than the default board that came with my system, but I had no idea just how much faster it would be, especially with iray scenes. When mentioning to Sean Kilbride, the technical marketing manager at NVIDIA for workstation product reviews, that this board would probably be great for games, he mentioned that this board is really a professional level board and would be overkill for games.

  I found that the Quadro 5000 graphics board significantly speeded up my workflow and made me much more productive, especially when I got to the test render phase of a project. I found I could quickly find where I needed to spend more time without having to do a complete render. The end result is that I can get more projects done in a shorter amount of time, which is great news when I got yet another deadline looming.

  You can find more information on 3ds Max and NVIDIA iray at the www.autodesk.com web site and more information on the Quadro line of professional graphics cards is available at www.nvidia.com/quadro.

MLAA: Efficiently Moving Antialiasing from the GPU to the CPU

Tue, 26 Jul 2011 12:54:47 +0000

Brought to you by the Intel® Visual Computing Developer Community

Download the Source Code »
Watch the Video »

Introduction

Efficient antialiasing techniques are an important tool of high-quality, real-time rendering. MSAA (Multisample Antialiasing) is the standard technique in use today, but comes with some serious disadvantages:

Incompatibility with deferred lighting, which is used more and more in real-time rendering;
High memory and processing cost, which makes its use prohibitive on some widely available platforms (such as the Sony Playstation* PS3* [Perthuis 2010]). This cost is also directly linked to the complexity of the scene rendered;
Inability to smooth non-geometric edges unless used in conjunction with alpha to coverage.

A new technique developed by Intel Labs called Morphological Antialiasing (MLAA) [Reshetov 2009] addresses these limitations. MLAA is an image-based, post-process filtering technique which identifies discontinuity patterns and blends colors in the neighborhood of these patterns to perform effective antialiasing. It is the precursor of a new generation of real-time antialiasing techniques that rival MSAA [Jimenez et al., 2011] [SIGGRAPH 2011].

This sample is based on the original, CPU-based MLAA implementation provided by Reshetov, with improvements to greatly increase performance. The improvements are:

Integration of a new, efficient, easy-to-use tasking system implemented on top of Intel® Threading Building Blocks (Intel® TBB).
Integration of a new, efficient, easy to use pipelining system for CPU onloading of graphics tasks.
Improvement of data access patterns through a new transposing pass.
Increased use of Intel® SSE instructions to optimize discontinuities detection and color blending.

The MLAA algorithm

In this section we present an overview of how the MLAA algorithm works; cf. [Reshetov 2009] for the fully detailed explanation. Conceptually, MLAA processes the image buffer in three steps:

Find discontinuities between pixels in a given image.
Identify U-shaped, Z-shaped, L-shaped patterns.
Blend colors in the neighborhood of these patterns.

The first step (find discontinuities) is implemented by comparing each pixel to its neighbors. Horizontal discontinuities are identified by comparing a pixel to its bottom neighbor, and vertical discontinuities by comparing with the neighbor on the right. In our implementation we compare color values, but any other approach that suits the application’s specificities is perfectly valid.

At the end of the first step, each pixel is marked with the horizontal discontinuity flag and/or the vertical discontinuity flag, if such discontinuities are detected. In the next step, we “walk” the marked pixels to identify discontinuity lines (sequences of consecutive pixels marked with the same discontinuity flag), and determine how they combine to form L-shaped patterns, as illustrated below:

Figure 1: MLAA processing of an image, with Z, U, and L-shapes shown on the original image on the left

The third and final step is to perform blending for each of the identified L-shaped patterns.

The general idea is to connect the middle point of the primary segment of the L-shape (horizontal green line in the figure below) to the middle point of the secondary segment (vertical green line—the connection line is in red). The connection line splits each pixel into two trapezoids; for each pixel the area of the corresponding trapezoid determines the blending weights. For example, in the figure below, we see that the trapezoid’s area for pixel c5 is 1/3; so the new color of c5 will be computed as 2/3 * (color of c5) + 1/3 * (color of c5’s bottom neighbor).

Figure 2: Computing blending weights

In practice, to ensure a smooth silhouette look, we need to minimize the color differences at the stitching positions of consecutive L-shapes. To achieve this, we slightly adjust the connection points on the L-shape segments around the middle position based on the colors of the pixels around the stitching point.

Sample usage

The camera can be moved around the scene using drag-and-click, and the mouse wheel can be used to zoom in or out. In addition, there are three blocks of controls on the right side of the sample’s window:

Figure 3: A screenshot of the sample in action

The first block controls the rendering of the scene:

Pause Scene toggles the scene’s animations on/off;
Show zoom box toggles the zoom box feature on/off, which allows the user to take a closer look at a pixel area to better compare the antialiasing techniques. The area of interest can be changed by right-clicking to the new area to examine;
Scene complexity simulates the effect of increasing scene complexity by using overdraw. The value (between 1 and 100, adjusted by the slider) indicates how many times the scene is rendered per frame (with overdraw).

The second block selects which antialiasing technique to apply: MLAA, MSAA (4x) or no antialiasing. This is of course to allow comparison of the techniques’ performance and quality (the default choice is “MLAA”);

The last block of controls is only available if the antialiasing technique used is MLAA, and controls how the algorithm should be run:

Copy/Map/Unmap only copies the color buffer from GPU memory to CPU memory and back, but won’t perform any MLAA processing between the two copy operations. This allows measurement of the impact of the copy operations on the overall performance of the entire algorithm;
CPU tasks pipelining turns the pipelining system on/off for CPU onloading of graphics tasks (the default is on) so that it is easy to see the benefit of pipelining;
Show found edges runs the first part of the algorithm (find discontinuities between pixels), but the blending passes are replaced by a debug pass, where a pixel is:
changed to solid green if a horizontal discontinuity has been found with its neighbor;
to solid blue if a vertical discontinuity has been found with its neighbor;
changed to solid red if both horizontal and vertical discontinuities have been found with its neighbors;
unchanged if no discontinuities have been found.

Figure 4. Gear tower scene with MLAA and the ”show found edges” option enabled

Sample Architecture

Without the pipelining optimizations, the sequence of events for each frame is:

ANIMATE AND RENDER TEST SCENE

MLAA STAGE (if MLAA is enabled)

Copy the color buffer where the scene was rendered to a staging buffer
Map the staging buffer for CPU-side access.
MLAA post-processing (the staging buffer is both input and output)
Unmap the staging buffer
Copy staging buffer back to the GPU-side color buffer
Render the zoombox (if applicable)

RENDER THE SAMPLE’s UI, PRESENT FRAME

Except for the “Perform the MLAA post-processing work” step (and again not considering the pipeline for now), all of these steps are implemented using standard Microsoft DirectX* methods.

The Tasking API

By nature the MLAA algorithm is easy to parallelize. For both the discontinuities detection and blending passes, we can process the color buffer in independent chunks (blocks of contiguous rows or columns). MLAA can be implemented using a task-based solution that automatically makes full use of all available CPU cores, while keeping the code core count agnostic

This sample uses a simple C-based tasking API that is implemented on top of the Intel® Threading Building Blocks (Intel® TBB) scheduler. This wrapper API was created to simplify the integration of the technique into existing codebases which already expose a similar tasking API (e.g. cross-platform game engines). An added benefit is the increased readability of the main source file MLAA.cpp.

The two important functions of the wrapper APIs are:

And:

The callback function has the following signature:

.MLAA requires a dependency graph of three consecutive tasksets:

The first taskset finds the discontinuities between pixels in the color buffer;
The second taskset performs the horizontal blending pass, and depends on the completion of the first taskset for the discontinuities information;
The third taskset performs the vertical blending pass, which depends on the completion of the first taskset for the discontinuities information, but also on the completion of the second taskset because of the transpose optimization (cf. corresponding section for details).

This dependency graph is expressed in MLAA.cpp as a sequence of three calls to CreateTaskSet; the taskset callback functions (implemented in MLAAPostProcess.cpp) being respectively MLAAFindDiscontinuitiesTask, MLAABlendHTask, MLAABlendVTask.

If pipelining is not enabled, we wait for the last taskset to complete by calling WaitForSet with the handle of the last taskset. When the call returns, the MLAA work for the frame is complete. Things are just a little bit more complicated when using the pipeline.

The CPU/GPU pipeline

To get the maximum performance from our implementation, we have to keep both CPU and GPU sides fully utilized. Due to the data dependencies between the tasksets, full utilization can be achieved by interleaving the processing of multiple frames, as shown in the diagram below:

Figure 5: Frames moving through the pipeline. The red blocks on the CPU main thread timeline illustrate that the main thread
will take MLAA work if the application gets CPU-bound.

In other words, we have to build a pipelining system, where a pipeline is a sequence of CPU stages (workloads) and GPU stages, and able to run multiple instances of said pipeline at the same time.

In our case, each instance of the pipeline corresponds to the processing of one frame. We have three steps:

Step 1 (GPU stage):

Animate and render test scene;
Copy color buffer to staging buffer (using asynchronous GPU-side CopyResource).

Step 2 (CPU stage):

Map the staging buffer;
Perform the MLAA post-processing work (staging buffer is both input and output).

Step 3 (GPU stage):

Unmap staging buffer and copy it back to the GPU-side color buffer;
Finish frame rendering (zoombox, UI) and present.

To implement this concept, we designed a simple Pipeline class. Each stage is represented by a PipelineFunction structure that specifies the function to be called, and the stage type. The callback function must have the following signature:

Where uInstance indicates which instance of the pipeline the function is being called from. A GPU stage uses a DirectX* query (of type D3D11_QUERY_EVENT) to signal its completion, but a CPU stage has to explicitly call the Pipeline’s class CompleteCPUWait method to signal completion:

Creating the pipeline instances requires a single call to the Init method. In our case, the code is:

Where g_NumPipelineInstances is the number of pipeline instances to create (3 in this case).

To run the pipelines, we call the Start method each frame:

Because the Start method returns the index of a completed instance of the pipeline, the third and last step doesn’t have to be called through the Pipeline class. The code of the last step executes right after the call to the Start method, indexing the data structures with the index returned by Start.

Pipelining does not rely on the tasking API WaitForSet call, since it is blocking and so does not allow pipelining to occur. The solution is to use a completion taskset—that is, a task that depends on the completion of all the MLAA tasks, and whose only work will be to call the CompleteCPUWait method.

Intel® SSE optimizations

The first pass of the MLAA algorithm identifies discontinuities between pixels. Conceptually, each pixel checks its color and compares it with its bottom neighbor (when looking for horizontal discontinuities) or right neighbor (when looking for vertical discontinuities) [Reshetov09, section 2.2.1].

In this sample, we have kept the simple discontinuity detection kernel of the reference implementation. A discontinuity exists between two pixels if at least one of their RGB color components differs by at least 16 (on the 0-255 scale of the RGBA8 format).

This definition works well in the sample and allows a very compact and efficient SIMD implementation. However, more complex approaches are possible, and sometimes necessary, to get the best possible results. For example, a pixel’s luminance could be used instead of directly comparing color values, and a variable threshold recomputed each frame to take into account the scene’s overall luminance and contrast [Luminance]. Depth values could be used to assist with edge detection or any kind of custom data to exclude specific zones of the color buffer from being processed.

Because this step is independent from the rest of the algorithm, and we are running on the CPU, any data from the program can be used directly as input data. The only limits to the detection algorithm are:

The fact that the algorithm must output a bit for each pixel indicating whether a discontinuity is detected;
The tradeoff between performance and quality/complexity;
The programmer’s imagination.

As with the original reference implementation [Reshetov09, section 2.4], we work with a RGBA8 render target format and the “vertical discontinuity” and “horizontal discontinuity” bit flags computed by MLAA are stored in-place in the two high bits of the pixel data (i.e., the two high bits of the alpha component which is unaffected by the MLAA blending operations). This keeps the implementation simple while helping with memory footprint, and allowing optimizations in the next steps of the algorithms (cf. “the transposing optimization” section below).

Because each pixel is 32 bits of data, and each pixel can be processed independently of the others in this step, we can process 4 pixels at a time using Intel® SSE intrinsics. The detection code is very short.

The reference implementation then proceeds to update the alpha component of each pixel, but used inefficient serial code to do so. We can optimize this sequence by using the following Intel® SSE code:

We also replaced the MixColors function (which computes a linear interpolation of two colors) to use a full Intel® SSE implementation.

The Transposing Optimization

The next task of the algorithm is to find “discontinuity lines”, i.e., sequences of consecutive pixels which are marked with the same discontinuity flag (horizontal flag for horizontal blending pass, vertical flag for vertical blending pass) while walking rows and columns of our color buffer [Reshetov09, section 2.1].

Because discontinuity lines tend to be short, intuition suggests (and profiling data confirms) that the most expensive part of this operation is scanning between discontinuity lines, i.e., scanning the large areas of consecutive pixels where no discontinuity flag is set.

The good news is that we can check 4 pixels at a time using the _mm_movemask_ps Intel® SSE intrinsic when the following conditions are met:

We are scanning 4 pixels stored at consecutive addresses;
The address of the first pixel is “Intel® SSE-aligned” (16 bytes alignment);
The discontinuity flag is stored as the high bit of the 32 bits of pixel data.

During the horizontal blending pass, (1) is true (we scan horizontal rows of pixels in the color buffer represented as a 2D linear array of pixel data); (2) is true almost all the time (remember that 16 bytes alignment is equivalent to “the index of the starting pixel in the buffer is a multiple of 4” as the buffer is properly aligned); and (3) is true as we chose bit 31 to represent the horizontal discontinuity flag.

If all conditions are true, we compute the flags:

And five outcomes are possible depending on the value of HFlags:

0 (the most common case by far): no discontinuity flag set in this group of 4 pixels, move to the next group of 4 pixels.
Bit 0 is set: discontinuity line starts at first pixel of the group.
Bit 1 is set: discontinuity line starts at second pixel of the group.
Bit 2 is set: discontinuity line starts at third pixel of the group.
Bit 3 is set: discontinuity line starts at fourth pixel of the group.

This optimization is part of the reference implementation and works well for the horizontal blending pass, but was impossible to apply for the vertical blending pass. As the code scans the buffer vertically, (1) is false (adjacent pixels in a column are not stored at consecutive addresses), and (3) is false as well (the vertical discontinuity flag is stored at bit 30 of the pixel data).

The problem with (3) is easy to work around by introducing a simple “shift left by one bit” operation if we are processing the vertical flags, transforming the code above to:

But (1) is still a problem. In addition, vertical scanning is extremely cache-unfriendly. Because of both of these issues, the vertical blending pass is about 3 times (300%!) as expensive as the horizontal blending pass in the reference implementation (as shown by profiling data).

Our solution to this issue is to make the vertical blending pass use the cache and Intel® SSE-friendly data access patterns of the horizontal blending pass by considering the color buffer as a matrix of pixels and transposing it between passes:

Perform horizontal blending pass
Transpose the (horizontally blended) color buffer
Perform vertical blending pass
Transpose back color buffer

This way the code for both blending passes is exactly the same (which adds the advantage of simplicity/readability), the only difference being which flag to scan for, and we benefit from all the optimizations and cache-friendliness of the horizontal pass.

In practice, the transpose operations are not implemented as separate passes/tasksets, but as the last part of their respective blending pass. This allows us to benefit from the cache “warmness”.

The two transpose operations are not free. The cost is the extra code execution time, one extra work buffer and a synchronization point between the horizontal and vertical passes (we must wait for all horizontal tasks to be done before we can start any of the vertical tasks, because we must wait for the color buffer to be fully transposed before starting the vertical pass work).

Even with these extra costs, the overall performance is significantly better than the previous approach. As expected, profiling data shows that both blending passes have equivalent performance.

Performance Results

MLAA performance was measured on the following two configurations:

Code name “Huron River” : Intel® Core™ i7-2820QM Processor (Intel® microarchitecture code name Sandy Bridge, 4 cores 8 threads @2.30 GHz) with GT2 processor graphics, 4 GB of RAM, Windows 7 Ultimate 64-bit Service Pack 1
Code name “Kings Creek”: Intel® Core™ i5-660 Processor (codename “Clarkdale”, 2 cores 4 threads @ 3.33 Ghz), with GMA HD Graphics (codename “Ironlake”), 2 GB of RAM, Windows 7 Ultimate 64-bit Service Pack 1

We measured the average frame rendering time of our sample as a function of scene complexity for the different antialiasing settings. We also measured the rendering time for the Copy/Map/Umap only mode to highlight the impact of the color buffer copy operations on the overall performance of the algorithm.

Results

The data shows that for the Huron River machine, the extra cost of MSAA 4x goes up linearly with the scene complexity (in fact, for all resolutions, the frame time when using MSAA 4x to render our test scene appears to be approximately the double of the frame time when no antialiasing is used). In contrast, the cost of MLAA appears more or less constant (around 4 ms/frame at 1280x800). This is consistent with our expectations, as unlike MSAA 4x, MLAA is executed only once per frame, regardless of scene complexity / number of draw calls.

We also observe that except for very low complexity values (i.e. less than ~5) MLAA always outperforms MSAA 4x, regardless of resolution, and that the difference in performance grows with complexity (because as noted above, the cost of MSAA 4x grows linearly with complexity when the cost of MLAA does not).

In the case of the Kings Creek machine, we can’t compare the cost of MLAA vs. MSAA 4x as the latter is not provided by the Ironlake hardware. The goal is then to determine if MLAA allows us to provide software antialiasing as an alternative with acceptable performance. Our measurements at the 1280x800 resolution show that again the cost of MLAA is largely independent of scene complexity, and the average value is ~7.5 ms (discarding the outlier data point at complexity = 100).

Interestingly, if we compare this result to a hypothetical MSAA 4x implementation with approximately the same performance profile than the Huron River one (frame time with MSAA 4x ~ 2x frame time with no antialiasing), we notice that again MLAA would outperform MSAA 4x for almost all complexity values (=4 in this case).

Table 1. Rendering times of our test scene on the King’s Creek machine.

Table 2. Rendering times of our test scene on the Huron River machine.

Figure 6: Rendering times of our test scene for each of the antialiasing techniques, as a function of scene complexity. The
bottom, flat curve is the difference between the “MLAA with pipeline on” curve, and the “No antialiasing” curve, measuring the
cost of using MLAA [Huron River, res. 1280x800]

Figure 7: Rendering times of our test scene for each of the antialiasing techniques, as a function of scene complexity. The
bottom, flat curve is the difference between the “MLAA with pipeline on” curve, and the “No anti aliasing” curve, measuring the
cost of using MLAA [Huron River, res. 1600x1200]

Figure 8: Rendering times of our test scene with MLAA on and off, as a function of scene complexity. The bottom, flat curve is
the difference between the “MLAA with pipeline on” curve, and the “No antialiasing” curve, measuring the cost of using MLAA
[Kings Creek, res. 1280x800]

References

[Reshetov 2009] RESHETOV, A. 2009. “Morphological Antialiasing”

[Perthuis 2010] PERTHUIS, C. 2010. MLAA in God of War 3. Sony Computer Entertainment America, PS3 Devcon, Santa Clara, July 2010.

[Jimenez et al., 2011] JIMENEZ, J., MASIA, B., ECHEVARRIA, J., NAVARRO, F. and GUTIERREZ, D. 2011. Practical Morphological Anti-Aliasing. In Wolfgang Engel, ed., GPU Pro 2. AK Peters Ltd.

[SIGGRAPH 2011] JIMENEZ, J., GUTIERREZ D., YANG, J., RESHETOV, A., DEMOREUILLE, P., BERGHOFF, T., PERTHUIS, C., YU, H., MCGUIRE, M., LOTTES, T., MALAN, H., PERSSON, E., ANDREEV, D. and SOUSA T. 2011. Filtering Approaches for Real-Time Anti-Aliasing. In ACM SIGGRAPH 2011 Courses.

[Luminance] Definition of luminance for CRT-like devices:

INTERNATIONAL COMMISSION ON ILLUMINATION. 1971. Recommendations on Uniform Color Spaces, Color Difference Equations, Psychometric Color Terms. Supplement No.2 to CIE publication No. 15 (E.-1.3.1), TC-1.3, 1971.

And for LCDs:

ITU-R Rec. BT.709-5. 2008. Page 18, items 1.3 and 1.4

Comparing Shadow Mapping Techniques with Shadow Explorer

Mon, 21 Mar 2011 17:26:13 +0000

This article brought to you by Intel® Visual Computing Developer Community

Code Sample Download Page

http://software.intel.com/en-us/articles/shadowexplorer/

Download Article

Download Comparing Shadow Mapping Techniques with Shadow Explorer [PDF 990KB]

Introduction

Most modern games use shadows to some degree. The challenge is to know which algorithm to use and what the tradeoffs are of different techniques-what combination of quality and performance is best suited for the game. This sample, Shadow Explorer, lets the user compare and contrast four different algorithms, adjust various parameters for each one, and observe the effects in real time. The shadow mapping algorithms presented are: simple, percentage closer filtered (PCF), variance (VSM), and exponential variance (EVSM).

Sample Usage

At the highest level, this sample allows the user to compare and contrast the quality and performance characteristics of the different shadowing algorithms. There are two scenes that each algorithm can be applied to: a city and a teapot. The city scene represents a "typical" setting that might be encountered in a game, consisting of a variety of object sizes, dimensions, and other characteristics. The teapot is a "worst case" scenario consisting primarily of thin primitives, casting shadows not only on themselves but also on curved surfaces and a flat plane.

Sample Architecture

Shadow Explorer is implemented as a Microsoft DirectX* application based on the DXUT framework. All four shadow mapping algorithms basically follow the same code path, except that different shaders are used and that VSM and EVSM run a filter over the shadow map before the scene is rendered. Listing 1 shows the basic flow of the source code for creating and using the shadow map. Every frame, the shadow map is created in RenderShadowMap() which builds a standard cascaded shadow map by running the SceneZ shader. As mentioned, for VSM and EVSM the shadow map is then filtered using the FilterV and FilterH shaders. After the shadow map is created, the scene is rendered with the shader SceneMain which, depending on the algorithm selected, will call a different version of IsNotInShadow(). This is done by simply including different .fxh files and recompiling the shaders when a new algorithm is selected. For example, if PCF is selected then IsNotInShadow() will call the version in filter_PCF.fxh.

static void RenderShadowMap(...)
{
   // For each cascade
   for( int i = 0; i < iLayers; ++i )
   {
	  // Set the Z only pass
	  SetRenderTargets( 0, NULL, g_d3d.pShadowMapDSV[i] );
	  g_d3d.pTechSceneZ->GetPassByIndex( 0 )->Apply( 0,
		 pd3dImmediateContext );
	  // Render the appropriate scene
	  g_pSelectedMesh->Render( pd3dImmediateContext );
   }
	
   // VSM and EVSM techniques require additional filters
   if(Filter_Type_VSM || Filter_Type_EVSM)
   {
	  // Use the shadow map just generated in the filters
	  g_d3d.pVarShadowTex->SetResource( g_d3d.pShadowMapSRV );

	  // Run the vertical filter
	  SetRenderTargets( 1, &pShadowMapRTV_vsm[0].p, NULL );
	  g_d3d.pTechFilter_V->GetPassByIndex( 0 )->Apply(...);
	  pd3dImmediateContext->Draw(...);

	  // Use the result of the vertical filter
	  g_d3d.pVarShadowTex->SetResource( g_d3d.pShadowMapSRV_vsm[0] );

	  // Run the horizontal filter
	  SetRenderTargets( 1, &pShadowMapRTV_vsm[1].p, NULL );
	  g_d3d.pTechFilter_H->GetPassByIndex( 0 )->Apply(...);
	  pd3dImmediateContext->Draw(...);
   }
}

static void RenderScene(...)
{
   RenderShadowMap(...)

   if(Filter_Type_VSM || Filter_Type_EVSM)
	  // Use the filtered shadow map
	  g_d3d.pVarShadowTex->SetResource( g_d3d.pShadowMapSRV_vsm[1] );
   else
	  g_d3d.pVarShadowTex->SetResource( g_d3d.pShadowMapSRV );
   g_pSelectedMesh->Render(...);
}

Listing 1 - Basic flow of application code in Shadow Explorer

The simplest technique presented is that of a basic, cascaded shadow map. This is obviously going to be the fastest method as it does the least amount of work, but leaves a lot to be desired in terms of quality. The next technique is PCF, which works off the same shadow map but samples the texture multiple times to alleviate aliasing problems and soften the shadow edges. Better results can be achieved by taking more samples, but with the expected tradeoff in performance. By using a non-uniformly distributed sample distribution, fewer samples need to be taken to achieve good results. Shadow Explorer uses a pre-calculated Halton sequence to calculate "random" offsets instead of a grid based offset. Even using the non-uniform method, PCF still requires a lot of sample look ups, making the algorithm fairly time consuming on the GPU.

The next technique, VSM, takes a different approach by storing both the depth and depth squared values into the shadow map. In Shadow Explorer, the shadow map is rendered normally and the depth squared value is added during a separate pass immediately after the z-pass is done. Additionally, this second pass runs a box filter over the data to soften the edges of the shadows. The main drawback to the VSM algorithm is the light bleeding effect when multiple occluders overlap each other and the ratio of their distances from the shadow receiver is high. This can be seen in Figure 1, where the shadow of the tall building in the background is outlined in light where it overlaps the shadows of the two smaller buildings in the foreground.

Figure 1 - VSM results in light bleeding artifacts

Depending on the scene, light bleeding may not be noticeable enough to be a problem. If it is, EVSM can be used to eliminate the problem, as can be seen in Figure 2. The main difference with VSM is that EVSM "warps" the depth value when storing and reading from the shadow map. This has the effect of bringing the relative distances of the occluders to the receiver closer together in order to minimize or eliminate light bleeding. Shadow Explorer does this in the shader function WarpDepth().

Figure 2 - EVSM fixes light bleeding

Conclusion

Shadow Explorer implements four shadow mapping techniques with runtime modification of various parameters so the user can compare both the performance and the quality of the methods. For each technique, care was taken to ensure optimal performance of the algorithms. Some of the optimizations done include using half floats instead of floats, unrolling loops in PCF, running the vertical filter before the horizontal filter, selecting the cascade Z interval instead of fitting to best cascade, and calculating shadows only for triangles oriented towards the light. Additional resources, including both binary and source code versions, are available for download at http://software.intel.com/en-us/articles/shadowexplorer.

Some areas for further development would be the addition of spot and point lights to demonstrate the flexibility of the shadow techniques. Animated objects could also be added to create an environment more representative of a game. Also, different algorithms could be used to more accurately place the split planes.

Controls

Shadow Explorer allows the user to change a variety of parameters to modify the shadow algorithms and change the shadow technique being used. Here is a list of the controls along with a brief description:

Toggle full screen - Turns full screen on/off
Change device - Change the d3d device
Scene drop down list - Change the scene being viewed
Align light to camera - Rotates the light along the light direction for optimal usage of shadow map space.
Algorithm drop down list - Change the filtering method being used
Shadow Map BPP - How many bits are used for each pixel in the shadow map
SM resolution - resolution of the shadow map
Filter size - The size of the filter used in PCF, VSM, and EVSM
Cascade layers - How many cascades are used
Cascade Factor - Adjusts how the split space is partitioned by split planes
Aperture - Spatial screen space filter aperture for PCF
Visualize Cascades - Visualize the different cascade levels
Visualize Light Space - Displays the scene from the light source's point of view.
Use Texture Array - Use texture array for cascades. Each cascade layer receives one texture, otherwise a texture atlas is used.
Z interval selection - How to determine which cascade a particular pixel is in. If enabled, then only the view space z distance is used. Otherwise, the world position of the pixel is tested against each cascade.
Deduce Z range - Calculates z-range from the current view. Otherwise uses the scene bounding box as the upper bound.
Deduce Res - Resolution of the render target used to calculate z range.
Downscale Factor - Every subsequent GPU deduce pass uses a render target texture of smaller size. This option determines the relation between texture sizes of subsequent passes.
Downscale Limit - Maximum texture size when Z-range deduction is performed on GPU. At some point it should be more efficient to pull the whole resource to the CPU and finish deduction there, saving several draw calls.

References

Williams, L. 1978. Casting curved shadows on curved surfaces. In Proc. SIGGRAPH, vol. 12, 270-274. http://portal.acm.or...n.cfm?id=807402
Donnelly, W. and Lauritzen, A. Variance shadow maps. In SI3D '06: Proceedings of the 2006 symposium on Interactive 3D graphics and games. 2006. pp. 161-165. New York, NY, USA: ACM Press. http://www.punkuser....m/vsm_paper.pdf
Lauritzen, Andrew and McCool, Michael. Layered variance shadow maps. Proceedings of graphics interface 2008, May 28-30, 2008, Windsor, Ontario, Canada. http://portal.acm.or...=GUIDE&dl=GUIDE
Isidoro, J. R. Shadow Mapping: GPU-based Tips and Techniques. Conference Session. GDC 2006. March 2006, San Jose, CA. http://developer.amd...adowMapping.pdf
The Halton Sequence. http://orion.math.ia...n_sequence.html
Engel, Wolfgang F. Section 4. Cascaded Shadow Maps. ShaderX 5, Advanced Rendering Techniques, Wolfgang F. Engel, Ed. Charles River Media, Boston, Massachusetts. 2006. pp. 197-206.

Read the Blog!

http://origin-software.intel.com/en-us/blogs/2011/02/23/another-gaming-and-graphics-sample-coming-to-a-download-near-you/

Real-Time Dynamic Fur on the GPU

Thu, 14 Oct 2010 14:30:22 +0000

Particle Systems
In computer graphics, particle systems represent an unstable side of nature: substances that donÆt have a surface and change rapidly over time, like rain, fire, fog, sparks. Every game engine has some kind of a sub-system for processing the particles, ideally exposing the full set of the functionality in the native particle editor. Until now, when OpenGL 3 and GPGPU tools gained wide support, particle engines were limited to 2¹⁰-2¹³ particles being processed simultaneously.  This article describes an advanced GPU-only particle engine structure. The technologies described are designed in terms of OpenGL core profile, but can be ported to Direct3D 10 as well. The reader should have a general understanding of particle systems, OpenGL-3 pipeline and modern high-level shading languages.

Traditional Approach
  struct Particle {    //attributes    Vector3 pos, speed;    float size; }; class ParticleSystem {    std::vector array;    void init(unsigned num);    void step(float delta); };  The step() which updates the particle state including death and birth events is called all the time. The particle data is sent to GPU memory (represented by the buffer object) either by the direct data transfer call (glBuffer[Sub]Data(), glTexImage*) or implicitly by releasing of the locked memory (glUnmapBuffer) on every frame.
Point Sprites
The GPU is designed to process triangles, so each particle is generally drawn as a triangle (3 vertices) or a quad (2 triangles, 4 vertices). Point Sprite extension provides the functionality to generate the quad automatically on GPU from a given center coordinate and a size (glPointSize). When drawing points (GL_POINTS) in GL3+ context, it is enabled by default. An obvious advantage here is a great reduction in bandwidth load: sending 1 point data instead of 3-4 per particle. A downside is: the quad is generated in screen space and so can not be rotated and projected by the GL pipeline. Some particle systems work with symmetric forms (smoke, snow, fluids) and fit the point sprite model pretty well. In order to draw more complex primitives (like a rotated triangle), one can simulate a virtual viewport inside the point sprite quad area, implementing the rotation/projection transformations in the fragment shader and discarding unused area pixels.

Modern Approach
Most particles donÆt require processing of self-interactions (rain, fire, fog). This makes the main particle processing loop iterations independent and thus easy to run in parallel. The GPU can handle such processing in a more efficient way compared to the CPU, as it was originally designed for parallel computations of vertices and fragments. One of the features that made its way into the GL 3.1 core profile is Transform Feedback. It allows one to turn off the graphics rasterizer (optionally) and write transformed vertex data back into buffer objects. By using this technology and matching the input/output format of the vertex shader one can process particle data entirely on GPU. It will be stored in buffer objects, available for drawing right away. The difference from the CPU approach here is the necessity in the separate output buffer, as itÆs not possible to transform the data in-place on the GPU. The sequence of operations goes as follows:

init A; A->B; B->A; ... draw A ... ; A->B; ...

It is even possible to initialize particle data state right on the GPU. There is a special shader value gl_VertexID that allows us to distinguish different particles in a shader (only the initial values are set, so no input data is given at this stage).

It is reasonable to do particle processing and drawing separately with different frequencies. Thus the rasterizer is disabled when updating the data (GL_RASTERIZER_DISCARD). However, there is an option to combine these stages by enabling the rasterizer - and letting the transformed attributes to be drawn on the spot. This way is a bit faster, but much less scalable for the general purpose. It can be used for demos.

Implementation Details

Birth and Death
It is impossible to allocate memory in the shader, so the chunk of memory which is initially allocated for the particle data storage should already contain all possible particles that can co-exist. The allocated number of particles = capacity of the system = maximum theoretical number of particles living at the same time. The boolean flag isAlive can be encoded into one of the particle attributes (like particle moment of birth, where value <0 means that the particle is dead). The decisions of birth and death are made right inside the vertex shader that is responsible for particle update processing:

if(isAlive)    //update() returns 0 for dead particles    isAlive = update(); else if( born_ready() )    //reset() gives initial state to a new-born    isAlive = reset();
Managing Multiple Emitters
If there are several emitters of the same particle format and capacity, they can share a single temporary buffer for data processing. It can be implemented by a ParticleManager class that contains the same data as each of the emitters E0,E1,E2,... The data flow is the following: Ei -> M (processing); Ei <> M (swap data)

Each emitter has the current copy of the particle data at any moment.

Conclusion
The proposed technology of particle systems processing provides a more efficient utilization of the GPU power, reduces the system bandwidth load and frees CPU memory for other tasks like physics and AI.
Simulating Fur Strands
In the previous page the technique of modeling particle systems on the GPU was covered. Now weÆll talk about an extension to that technique that allows modeling of fur strands. In addition to the previous pre-requisites, we expect the readerÆs understanding of Texture Buffer Objects, Render Target concept and spatial transformations, including quaternion operations. The Blender representation of the strand that needs to be simulated consists of the fixed number of joints, connected by straight lines. A brief description can be found in Particle System and Fur Rendering, but the best source of knowledge is Blender itself. This model was chosen in order to have a possibility to edit hair in Blender and then easily move it into an engine. The export script into the KRI engine format has support for particles and hair, including physics and rendering parameters.

Definition: fur layer - a set of points joint (see picture) across all of the strands of a given surface.

Emit Sources
General particle data source is a mesh. Its geometry can be accessed directly from the shader that initializes the particle data on the birth event. This fact allows efficient emission of particles from the vertices or even the whole surface of a mesh, which can also be skinned and morphed at the same time.
Emit From Vertices
The way to the vertices lies through the Texture Buffer Object concept. It is required to bind the mesh data as TBO and sample from it in the particle processing shader, extracting the position and orientation (see Appendix A)  of a randomly chosen vertex. Once the vertex is extracted, it is possible to apply the initial direction of a particle in the constructed basis to get the speed and copy the position.  uniform samplerBuffer unit_vertex, unit_quat; int cord = int( random() * textureSize(unit_vertex) ); vec3 pos  = texelFetch(unit_vertex, cord).xyz; vec4 quat = texelFetch(unit_quat, cord);
Emit From Faces
In order to access the mesh surface data in the random pattern, it is necessary to bake it first into a UV texture. For each texel in the UV space the position and orientation (as a quaternion) of the point in the world space are stored, writing into 2 Rgba8 textures. In the vertex shader, the position is constructed from the texture coordinate:  in vec2 at_tex0;  //texture coordinate gl_Position = vec4( 2.0*at_tex0 - vec2(1.0), 0.0,1.0 );  When the textures are ready, a random sample is taken from them in the same way it was done for the vertices emitter case. Note: UV texture defines the amount and density of the grown fur, hence giving full control to the artist.

Fur Extension
Each fur layer is simulated as a separate particle emitter, all layers belong to the same particle manager. Capacity of that system equals the number of fur strands that need to be modelled. The strand joint particle is required to store at least position and speed attributes.  There are 3 issues arising from this point of view:

Making a Layer Dependent on Previous Layers
Definition: fake attribute - a particle attribute whose value is not stored per particle instance but rather taken from some external source (e.g. another particle emitter). The previous layer position is required in order to determine the current strand direction.

One layer before previous position is required for estimating the æstraightÆ direction and thus the current deviation. Both Pos(L-1) and Pos(L-2) are transmitted as fake attributes for processing of the current layer, where:

Pos(x) = position of the strand joint of layer æ[i]xÆ,
æ[i]Læ = current layer id (0

A Super Simple Method for Creating Infinite Scenery

Wed, 15 Sep 2010 11:53:45 +0000
My favorite kinds of games to play have always been realistic 3D vehicle simulations such as flight and ship simulations. So when I began game programming several years ago, I decided to come up with a way to create infinite scenery that was quick and easy to implement and did not require a lot of static geometry such as terrain data. Having such a method makes it easy to test camera and physics code because it does not require a terrain or water engine of any kind and you do not need to be worried about the size of the virtual world you will be testing in. Of course, the most common method is to simply create a very large rectangle to represent the ground and continually wrap the position of the vehicle to the opposite side every time it passes over one of the sides of the rectangle. This is infinite on the horizontal plane but when going vertically, eventually, the sides of the rectangle become visible because the rays of the camera frustum become long enough to go past the edges of the rectangle. So I came up with another method that was infinite in both the horizontal and vertical planes and requires no static geometry at all. In addition, it is very easy to implement and requires very little code and only simple math. The idea is to use the footprint of the camera, or in other words, the intersection of the camera’s frustum with the ground plane.

The simplest way to do this is to take our static ground rectangle and make it dynamic, resizing it as the camera moves or gets higher in altitude. But this means the rectangle will be much larger than it really needs to be most of the time. In this article, we will look at a more precise method that calculates the exact shape of the camera’s intersection with a flat ground plane preventing the need for view frustum culling or screen clipping.

Finding the camera’s footprint
In order to use this method, we need to have a camera position and orientation, as well as the dimensions of the viewport and the field-of-view angle.

The first step in the process is to calculate a rectangle that sits in front of the camera, is oriented according to the current orientation of the camera, and has the same width and height as the viewport.

Fig. 1 demonstrates this rectangle.

In order to calculate the corner points of the rectangle, we use the “side” and “up” unit vectors of the camera. In the demo, a camera class is used that stores these vectors but they can also be pulled out of the camera matrix if that is how the camera is represented instead. To get the y coordinate of each corner point we multiply the “up” vector by half the height of the viewport and to get the x coordinate we multiply the “side” vector by half the width of the viewport. In order to calculate the z coordinate we must know the “view distance”, or the distance along the “lookat” vector at which the frustum has the same dimensions as the viewport. Looking again at Fig. 1, we see that the rectangle that we are trying to create represents a single slice of the camera frustum. This slice can be made at any distance from the camera eye point, but at only one distance will the slice have the same width and height as the viewport. This distance is the “view distance”.

In order to calculate the “view distance” we use a simple formula based on the tangent of the field-of-view angle. In Fig.2 we have defined a triangle, shown in green, where the base, b, has the same length as the “view distance”. The height, h, is equal to half the viewport width and is parallel to the “side” vector. Since b bisects the field-of-view angle, then theta is equal to half the field-of-view angle.

Since we know that:

tan( theta ) = h / length of b

then if we solve for length of b we have

length of b = h / tan( theta )

So “view distance” is equal to the length of b, or half viewport width / half fov angle.

In the demo source code, the “drawHorizon” method of the RendererBase class contains the code for implementing these calculations:

  void RendererBase::drawHorizon( QCamera &qcamera, const ViewPort &viewport ) {    D3DXVECTOR3 lookAt = qcamera.lookAt();    D3DXVECTOR3 side = qcamera.side();    D3DXVECTOR3 up = qcamera.up();    D3DXVECTOR3 eye = qcamera.eye();    std::vector clippedPts, pts;    pts.reserve(4);    D3DXVECTOR3 viewC = lookAt * viewport.viewDistance();    float scrHalfW = viewport.screenWidth() * 0.5f;    float scrHalfH = viewport.screenHeight() * 0.5f;    // UL    pts.push_back((up * scrHalfH) + side * (-scrHalfW) + viewC); // UR    pts.push_back((up * scrHalfH) + (side * scrHalfW) + viewC);    //LR      pts.push_back((up * -(scrHalfH)) + (side * scrHalfW) + viewC);    //LL      pts.push_back((up * -(scrHalfH)) + (side * (-scrHalfW)) + viewC);  So now we will have four points contained in a vector and which represent the corners of our rectangle. The next step is to clip this rectangle below a horizontal plane that is parallel to the ground plane (i.e. all points have the same y-coord) and is slightly below the camera eye point. Why does it need to be below the camera eye and not exactly at the same height? This will be explained later but for now let us look at the clipping process. We will use the Cohen-Sutherland algorithm [LaMothe2003] , but pretty much any clipping algorithm will do. The basic idea is to look at each edge in the rectangle one-by-one and check to see if one of three cases is true: both end points are above the plane, both are below the plane, or one is above and the other is below. If both are above, we remove the edge from the list. If both are below we keep it in the list. If one is above, then we move this point down to the clipping plane and keep the other point the same. The rest of the “drawHorizon” method shows how this is done in code:  int v1,v2, cp = 0; float x1,y1,z1,x2,y2,z2,newx,newz,m; float ymax = -1.0f; v1=3; for (v2=0; v2 < 4; v2++) {    x1 = pts[v1].x;    y1 = pts[v1].y;    z1 = pts[v1].z;    x2 = pts[v2].x;    y2 = pts[v2].y;    z2 = pts[v2].z;    if ((y1 <= ymax) && (y2 <= ymax))    {   clippedPts.push_back( pts[v2] );    }    else if ((y1 > ymax) && (y2 > ymax))    {   // completely above    }    if ((y1 <= ymax) && (y2 > ymax))    {   if (x1 != x2)           {                              m= (y2-y1) / (x2-x1);                              newx = x1 + ((ymax - y1) / m);   }   else                      newx = x1;               if (z1 != z2)   {          m= (y2-y1) / (z2-z1);                              newz = z1 + ((ymax - y1) / m);   }   else          newz = z1;   clippedPts.push_back( D3DXVECTOR3(newx,ymax,newz) );    }    if ((y1 > ymax) && (y2 <= ymax)) {   if (x1 != x2)           {                              m= (y2-y1) / (x2-x1);                              newx = x1 + ((ymax - y1) / m);   }   else                      newx = x1;               if (z1 != z2)   {          m= (y2-y1) / (z2-z1);                              newz = z1 + ((ymax - y1) / m);   }   else          newz = z1;   clippedPts.push_back( D3DXVECTOR3(newx,ymax,newz) );   clippedPts.push_back( pts[v2] );    }    v1=v2; } cp = clippedPts.size(); if (cp == 0) {    return; }  We now have a vector containing a set of points which represent our clipped rectangle. This vector can contain anywhere from zero to five points depending on how the camera is oriented. Fig. 3 shows a diagram of this rectangle. The green area is our clipped rectangle, which represents the horizon image. The blue horizontal line represents our clipping plane. You can now see that if we draw lines from the camera eye to each of the rectangle corners, and then extend these lines until they intersect the ground plane, we will end up with a set of points which represent the camera’s footprint on the ground. This is why we needed to make the clipping plane slightly below our camera, because if it were at the same height, the top two lines would be parallel to the ground and would extend forever, never intersecting the ground. In the clipping code above, “yMax” represents the distance between the camera eye point and the clipping plane. In the demo, this is set to -1.0 but it can be any number less than zero. It is just a matter of visual preference, the further below the camera the clipping plane is, the further down in the viewport the horizon image will appear. Another option would be to define an epsilon value in our clipping code that would represent the minimum y-distance any point can be from the clipping plane. Then we would no longer need to move the clipping plane below the camera.

The code for extending the edges and converting into world space is shown below:

  for (unsigned int i=0; i < cp; i++) {    D3DXVec3Normalize( &footPrintVertex, &clippedPts[i] );    scale = fabs(footPrintVertex.y == 0 ? 0 : eye.y / footPrintVertex.y);    footPrintVertex *= scale;    footPrintVertex.x += eye.x;    footPrintVertex.y += eye.y;    footPrintVertex.z += eye.z;    ….    …. }  We simply go through each point in our vector of clipped points and normalize them. Then we scale each one by the height of the camera. Then we add the position of the camera eye and we now have our camera footprint in world space coordinates! Now all we have to do is render them as a triangle fan.
Rendering the camera’s footprint
In order to render our camera footprint, we need to consider what kind of projection we are using in our transformation pipeline.

If we are using a perspective projection, it is quite simple to render the footprint since it is already in world space coordinates. We simply pass these points into the API as a triangle fan and voila, we now have a horizon on screen. Of course, we should probably add textures to the horizon so we can have some ground detail which will give us a sense of altitude. This is also quite easy since all we have to do is scale the x and z coordinates of our points and we will get our tiled texture coordinates. This is shown in the “renderCameraFootprint” method of the PerspectiveRenderer class, included in the demo:

    geometryVertices[i].position = footPrintVertex; geometryVertices[i].tu = footPrintVertex.x / 200.0f; geometryVertices[i].tv = footPrintVertex.z / 200.0f; geometryVertices[i].ts = footPrintVertex.x / 15000.0f; geometryVertices[i].tt = footPrintVertex.z / 15000.0f;  The perspective render version can be seen in the viPerspective.exe application. We will use two textures to render the ground plane. This is because the higher you go, the smaller the texture tiles become until they are almost not visible. The second texture will be blended with the first and will have a larger tile size so it will become more visible as we go higher, providing new detail at higher altitudes. Now, if we are working in an environment where we will not be using any kind of a projection matrix, we will have to do a little more work since our footprint is in three-dimensional world coordinates and we need it to be in two-dimensional screen coordinates. Fortunately, the calculations are pretty easy, we just add perspective to our coordinates by dividing x and y by z, and then we scale by the view distance to convert them to viewport space. Since we are using a right-handed coordinate system in our demo, we use the negative view distance so that when the camera is pointing in the negative z direction the x coordinate will go from negative to positive when going from the left side of the viewport to the right side. Here is the code from the “renderCameraFootPrint” method of the OrthographicRenderer class from the demo:

  geometryVertices[i].x = ( (geometryVertices[i].x / geometryVertices[i].z) * -viewDist ); geometryVertices[i].y = ( (geometryVertices[i].y / geometryVertices[i].z) * -viewDist ); geometryVertices[i].z = 0.0f;  The main drawback is that we must also process our texture coordinates further. This will vary according to the platform and API that you are programming on. On the Playstation 2 it is relatively easy, we just calculate homogeneous texture coordinates s,t,q, like so:
s = geometryVertices[I].tu / footPrintVertex.z;
t = geometryVertices[I].tv / footPrintVertex.z;
q = 1.0f / footPrintVertex.z;
We then specify that perspective-correct texturing is enabled when rendering. I used this technique to create a skyplane cloud layer for a Playstation 2 flight simulator demo. Here is a screenshot:

YouTube Video Link Since the purpose of this article is to demonstrate the theory behind the method instead of any particular implementation details for specific graphics API’s, we will not go into how to do this on other platforms besides the PS2. In Direct3D and OpenGL it is probably best to stick to perspective projection since they will calculate the projection matrix for you and handle the texture coordinates as is.

But even without texturing, this technique can still be useful for creating an artificial horizon, such as on the instrument panel on airplanes, or for rendering the HUD of jet fighters. It is demonstrated in the viOrthographic.exe demo which uses an orthographic projection instead of a perspective projection since we do not need our perspective to be calculated by the API.

Adding Terrain
So great, we can now create a pretty nice looking horizon no matter where the camera is located or how it is oriented. However, so far we have only been able to render completely flat ground. What if we want to have some hills and mountains? Well, since we know what our footprint is on the two-dimensional ground plane, if we divide up the plane into a regular grid of equally sized rectangular “cells”, then all we have to do is figure out which cells overlap our footprint polygon and draw those cells using a heightmap to create three dimensional triangles and thus our terrain.

So how do we figure out which cells are overlapping? Well, if we think of our cells as pixels on a two-dimensional drawing surface, and our footprint polygon as a graphics primitive, then we can see how the scan conversion algorithm can be used, since it is the most common way to render graphics primitives onto a drawing surface.

The idea is that we create scan lines that have the same size as our grid cells and we scan along the line from one edge of the polygon that intersects the scan line to the next edge that intersects the scan line. All grid cells on the same scan line that lie between the two edges will be found by the scan.

If we look at Fig. 4, we can see how this technique works. On the first scan we start at cells that have a z-coordinate of -4. We then scan from left to right and find that only one grid cell intersects, cell (1, -4). We then increment our scan line to z-coordinate –3, and scan again starting at edge A and going to edge B. This time we find grid cells (0,-3) to (2,-3). We continue the process until we have scanned all of the edges. You can see how this is done in the demo source in the files ScanConvert.h and ScanConverter.cpp. There are many books and resources that explain scan conversion, especially ones that deal with software rasterization where scan conversion is quite often used.

The one drawback to this technique is that as the camera gets higher and higher, the footprint gets very large and thus requires a lot of grid cells to cover it. Eventually, the frame rate will drop off by quite a bit. In the demo, this problem is alleviated by increasing the size of the grid cells as the camera goes up. So when it reaches 20 units in altitude, the cell size is doubled to 2000 units, and when it reaches 40, it is quadrupled to 4000, and so on. The heightmap tile size however remains the same so that the terrain does not get wider, it just gets courser with less detail.

The viTerrain.exe demo demonstrates this technique in action.

As you can see there is quite a lot of popping so this method of optimization is only useful for demonstration purposes but it does help to improve the performance of the demo so that we can really get a good idea of how the scan conversion works. If you really want to get the best performance when using this technique with terrain, you should certainly use one of the common terrain LOD techniques to render each grid cell, such as ROAM, geomorphing, geomipmapping, etc. That is left as an exercise for the reader.

Having the camera’s world-space footprint can also be helpful when you do not need infinite terrain, but instead have a fixed set of geometry to represent your scenery. An example would be when you have several quadtrees to render terrain at different locations in world space. The camera footprint can be used in view-culling to determine if any of the geometry is visible from the camera. You simply need to perform 2D intersection tests between the footprint polygon and the geometry bounding boxes.

Conclusion
So there you have it, infinite scenery in both the horizontal and vertical dimensions with an extremely precise camera footprint that can be easily rendered in a very efficient manner.

Happy Coding,
Gabriel T. Delarosa
July 3, 2010

Bibliography:
Lamothe, Andre, “Tricks of the 3D Game Programming Gurus“ (Sams Publishing, 2003, ISBN: 0-672-31835-0 )

A Simple and Practical Approach to SSAO

Tue, 25 May 2010 13:41:25 +0000
Global illumination (GI) is a term used in computer graphics to refer to all lighting phenomena caused by interaction between surfaces (light rebounding off them, refracting, or getting blocked), for example: color bleeding, caustics, and shadows. Many times the term GI is used to refer only to color bleeding and realistic ambient lighting. Direct illumination – light that comes directly from a light source – is easily computed in real-time with today´s hardware, but we can´t say the same about GI because we need to gather information about nearby surfaces for every surface in the scene and the complexity of this quickly gets out of control. However, there are some approximations to GI that are easier to manage. When light travels through a scene, rebounding off surfaces, there are some places that have a smaller chance of getting hit with light: corners, tight gaps between objects, creases, etc. This results in those areas being darker than their surroundings.

This effect is called ambient occlusion (AO), and the usual method to simulate this darkening of certain areas of the scene involves testing, for each surface, how much it is “occluded” or “blocked from light” by other surfaces. Calculating this is faster than trying to account for all global lighting effects, but most existing AO algorithms still can’t run in real-time.

Real-time AO was out of the reach until Screen Space Ambient Occlusion (SSAO) appeared. SSAO is a method to approximate ambient occlusion in screen space. It was first used in games  by Crytek, in their “Crysis” franchise and has been used in many other games since. In this article I will explain a simple and concise SSAO method that achieves better quality than the traditional implementation.

The SSAO in Crysis

Prerequisites

The original implementation by Crytek had a depth buffer as input and worked roughly like this: for each pixel in the depth buffer, sample a few points in 3D around it, project them back to screen space and compare the depth of the sample and the depth at that position in the depth buffer to determine if the sample is in front (no occlusion) or behind a surface (it hits an occluding object). An occlusion buffer is generated by averaging the distances of occluded samples to the depth buffer. However this approach has some problems (such as self occlusion, haloing) that I will illustrate later.

The algorithm I describe here does all calculations in 2D, no projection is needed. It uses per-pixel position and normal buffers, so if you´re using a deferred renderer you have half of the work done already. If you´re not, you can try to reconstruct position from depth or you can store per-pixel position directly in a floating point buffer. I recommend the later if this is your first time implementing SSAO as I will not discuss position reconstruction from depth here. Either way, for the rest of the article I´ll assume you have both buffers available. Positions and normals need to be in view space.

What we are going to do in this article is exactly this: take the position and normal buffer, and generate a one-component-per-pixel occlusion buffer. How to use this occlusion information is up to you; the usual way is to subtract it from the ambient lighting in your scene, but you can also use it in more convoluted or strange ways for NPR (non-photorealistic) rendering if you wish.

Algorithm

Given any pixel in the scene, it is possible to calculate its ambient occlusion by treating all neighboring pixels as small spheres, and adding together their contributions. To simplify things, we will work with points instead of spheres: occluders will be just points with no orientation and the occludee (the pixel which receives occlusion) will be a  pair. Then, the occlusion contribution of each occluder depends on two factors:
Distance “d” to the occludee.
Angle between the occludee´s normal “N” and the vector between occluder and occludee “V”.
With these two factors in mind, a simple formula to calculate occlusion is: Occlusion = max( 0.0, dot( N, V) ) * ( 1.0 / ( 1.0 + d ) )

The first term, max( 0.0, dot( N,V ) ), works based on the intuitive idea that points directly above the occludee contribute more than points near it but not quite right on top. The purpose of the second term ( 1.0 / ( 1.0 + d ) ) is to attenuate the effect linearly with distance. You could choose to use quadratic attenuation or any other function, it´s just a matter of taste.


The algorithm is very easy: sample a few neighbors around the current pixel and accumulate their occlusion contribution using the formula above. To gather occlusion, I use 4 samples (<1,0>,<-1,0>,<0,1>,<0,-1>) rotated at 45º and 90º, and reflected using a random normal texture.

Some tricks can be applied to accelerate the calculations: you can use half-sized position and normal buffers, or you can also apply a bilateral blur to the resulting SSAO buffer to hide sampling artifacts if you wish. Note that these two techniques can be applied to any SSAO algorithm.

This is the HLSL pixel shader code for the effect that has to be applied to a full screen quad:

sampler g_buffer_norm; sampler g_buffer_pos; sampler g_random; float random_size; float g_sample_rad; float g_intensity; float g_scale; float g_bias; struct PS_INPUT { float2 uv : TEXCOORD0; }; struct PS_OUTPUT { float4 color : COLOR0; }; float3 getPosition(in float2 uv) { return tex2D(g_buffer_pos,uv).xyz; } float3 getNormal(in float2 uv) { return normalize(tex2D(g_buffer_norm, uv).xyz * 2.0f - 1.0f); } float2 getRandom(in float2 uv) { return normalize(tex2D(g_random, g_screen_size * uv / random_size).xy * 2.0f - 1.0f); } float doAmbientOcclusion(in float2 tcoord,in float2 uv, in float3 p, in float3 cnorm) { float3 diff = getPosition(tcoord + uv) - p; const float3 v = normalize(diff); const float d = length(diff)*g_scale; return max(0.0,dot(cnorm,v)-g_bias)*(1.0/(1.0+d))*g_intensity; } PS_OUTPUT main(PS_INPUT i) { PS_OUTPUT o = (PS_OUTPUT)0; o.color.rgb = 1.0f; const float2 vec[4] = {float2(1,0),float2(-1,0), float2(0,1),float2(0,-1)}; float3 p = getPosition(i.uv); float3 n = getNormal(i.uv); float2 rand = getRandom(i.uv); float ao = 0.0f; float rad = g_sample_rad/p.z; //**SSAO Calculation**// int iterations = 4; for (int j = 0; j < iterations; ++j) {   float2 coord1 = reflect(vec[j],rand)*rad;   float2 coord2 = float2(coord1.x*0.707 - coord1.y*0.707,   coord1.x*0.707 + coord1.y*0.707);      ao += doAmbientOcclusion(i.uv,coord1*0.25, p, n);   ao += doAmbientOcclusion(i.uv,coord2*0.5, p, n);   ao += doAmbientOcclusion(i.uv,coord1*0.75, p, n);   ao += doAmbientOcclusion(i.uv,coord2, p, n); } ao/=(float)iterations*4.0; //**END**// //Do stuff here with your occlusion value Ã¢aoÃ¢: modulate ambient lighting, write it to a buffer for later //use, etc. return o; }

The concept is very similar to the image space approach presented in “Hardware Accelerated Ambient Occlusion Techniques on GPUs” [1] the main differences being the sampling pattern and the AO function. It can also be understood as an image-space version of “Dynamic Ambient Occlusion and Indirect Lighting” [2] Some details worth mentioning about the code:
The radius is divided by p.z, to scale it depending on the distance to the camera. If you bypass this division, all pixels on screen will use the same sampling radius, and the output will lose the perspective illusion.
During the for loop, coord1 are the original sampling coordinates, at 90º. coord2 are the same coordinates, rotated 45º.
The random texture contains randomized normal vectors, so it is your average normal map. This is the random normal texture I use:



It is tiled across the screen and then sampled for each pixel, using these texture coordinates:

g_screen_size * uv / random_size

Where “g_screen_size” contains the width and height of the screen in pixels and “random_size” is the size of the random texture (the one I use is 64x64). The normal you obtain by sampling the texture is then used to reflect the sampling vector inside the for loop, thus getting a different sampling pattern for each pixel on the screen. (check out “interleaved sampling” in the references section)
At the end, the shader reduces to iterating through some occluders, invoking our AO function for each of them and accumulating the results. There are four artist variables in it:
g_scale: scales distance between occluders and occludee.
g_bias: controls the width of the occlusion cone considered by the occludee.
g_sample_rad: the sampling radius.
g_intensity: the ao intensity.
Once you tweak the values a bit and see how the AO reacts to them, it becomes very intuitive to achieve the effect you want.

Results


a) raw output, 1 pass 16 samples b] raw output, 1 pass 8 samples c) directional light only d) directional light – ao, 2 passes 16 samples each.

As you can see, the code is short and simple, and the results show no self occlusion and very little to no haloing. These are the two main problems of all the SSAO algorithms that use only the depth buffer as input, you can see them in these images:

♦
The self-occlusion appears because the traditional algorithm samples inside a sphere around each pixel, so in non-occluded planar surfaces at least half of the samples are marked as ‘occluded’. This yields a grayish color to the overall occlusion. Haloing causes soft white edges around objects, because in these areas self-occlusion does not take place. So getting rid of self-occlusion actually helps a lot hiding the halos.

The resulting occlusion from this method is also surprisingly consistent when moving the camera around. If you go for quality instead of speed, it is possible to use two or more passes of the algorithm (duplicate the for loop in the code) with different radiuses, one for capturing more global AO and other to bring out small crevices. With lighting and/or textures applied, the sampling artifacts are less apparent and because of this, usually you should not need an extra blurring pass.

Taking it further

I have described a down-to-earth, simple SSAO implementation that suits games very well. However, it is easy to extend it to take into account hidden surfaces that face away from the camera, obtaining better quality. Usually this would require three buffers: two position/depth buffers (front/back faces) and one normal buffer.  But you can do it with only two buffers: store depth of front faces and back faces in red and green channels of a buffer respectively, then reconstruct position from each one. This way you have one buffer for positions and a second buffer for normal.

These are the results when taking 16 samples for each position buffer:

left: front faces occlusion, right: back faces occlusion

To implement it just and extra calls to “doAmbientOcclusion()” inside the sampling loop that sample the back faces position buffer when searching for occluders. As you can see, the back faces contribute very little and they require doubling the number of samples, almost doubling the render time. You could of course take fewer samples for back faces, but it is still not very practical.

This is the extra code that needs to be added:

inside the for loop, add these calls:

ao += doAmbientOcclusionBack(i.uv,coord1*(0.25+0.125), p, n); ao += doAmbientOcclusionBack(i.uv,coord2*(0.5+0.125), p, n); ao += doAmbientOcclusionBack(i.uv,coord1*(0.75+0.125), p, n); ao += doAmbientOcclusionBack(i.uv,coord2*1.125, p, n);

Add these two functions to the shader:

float3 getPositionBack(in float2 uv) { return tex2D(g_buffer_posb,uv).xyz; } float doAmbientOcclusionBack(in float2 tcoord,in float2 uv, in float3 p, in float3 cnorm) { float3 diff = getPositionBack(tcoord + uv) - p; const float3 v = normalize(diff); const float d = length(diff)*g_scale; return max(0.0,dot(cnorm,v)-g_bias)*(1.0/(1.0+d)); }

Add a sampler named “g_buffer_posb” containing the position of back  faces. (draw the scene with front face culling enabled to generate it)  Another small change that can be made, this time to improve speed  instead of quality, is adding a simple LOD (level of detail) system to  our shader. Change the fixed amount of iterations with this:

int iterations = lerp(6.0,2.0,p.z/g_far_clip);

The variable  “g_far_clip” is the distance of the far clipping plane, which must be passed to the shader. Now the amount of iterations applied to each pixel depends on distance to the camera. Thus, distant pixels perform a coarser sampling, improving performance with no noticeable quality loss. I´ve not used this in the performance measurements (below), however.

Conclusion and Performance Measurements

As I said at the beginning of the article, this method is very well suited for games using deferred lighting pipelines because it requires two buffers that are usually already available. It is straightforward to implement, and the quality is very good. It solves the self-occlusion issue and reduces haloing, but apart from that it has the same limitations as other screen-space ambient occlusion techniques:  Disadvantages:
Does not take into account hidden geometry (especially geometry outside the frustum).
The performance is very dependent on sampling radius and distance to the camera, since objects near the front plane of the frustum will use bigger radiuses than those far away.
The output is noisy.
Speed wise, it is roughly equal to a 4x4 Gaussian blur for a 16 sample implementation, since it samples only 1 texture per sample and the AO function is really simple, but in practice it is a bit slower. Here´s a table showing the measured speed in a scene with the Hebe model at 900x650 with no blur applied on a Nvidia 8800GT:

Settings FPS SSAO time (ms)
High (32 samples front/back) 150 3.3
Medium (16 samples front) 290 0.27
Low (8 samples front) 310 0.08

In these last screenshots you can see how this algorithm looks when  applied to different models. At highest quality (32 samples front and  back faces, very big radius, 3x3 bilateral blur):

At lowest quality (8 samples front faces only, no blur, small radius):

It is also useful to consider how this technique compares to ray-traced AO. The purpose of this comparison is to see if the method would converge to real AO when using enough samples.

Left: the SSAO presented here, 48 samples per pixel (32 for front faces and 16 for back faces), no blur. Right: Ray traced AO in Mental Ray. 32 samples, spread = 2.0, maxdistance = 1.0; falloff = 1.0.

One last word of advice: don´t expect to plug the shader into your pipeline and get a realistic look automatically. Despite this implementation having a good performance/quality ratio, SSAO is a time consuming effect and you should tweak it carefully to suit your needs and obtain the best performance possible. Add or remove samples, add a bilateral blur on top, change intensity, etc. You should also consider if SSAO is the way to go for you. Unless you have lots of dynamic objects in your scene, you should not need SSAO at all; maybe light maps are enough for your purpose as they can provide better quality for static scenes.

I hope you will benefit in some way from this method. All code included in this article is made available under the MIT license

References

[1] Hardware Accelerated Ambient Occlusion Techniques on GPUs
(Perumaal Shanmugam)  [2] Dynamic Ambient Occlusion and Indirect Lighting
(Michael Bunnell)

[3] Image-Based Proxy Accumulation for Real-Time Soft Global Illumination
(Peter-Pike Sloan, Naga K. Govindaraju, Derek Nowrouzezahrai, John Snyder)

[4] Interleaved Sampling
(Alexander Keller, Wolfgang Heidrich)

Crytek´s Sponza rendered at 1024x768, 175 fps with a directional light.

The same scene rendered at 1024x768, 110 fps using SSAO medium settings: 16 samples, front faces, no blur. Ambient lighting has been multiplied by (1.0-AO).

The Sponza model was downloaded from Crytek's website.

Deferred Rendering Demystified

Mon, 01 Mar 2010 16:46:32 +0000
This article is a design article about implementing deferred rendering. The motive behind it is that while there have been many articles and presentations about the concepts behind deferred rendering (for example, the article about deferred rendering in Killzone 2), there is very little information about how to approach it from a design standpoint. This article aims to do just that.

The article is accompanied by code that implements a deferred rendering framework, and is somewhat a development journal of that framework. You can get the code from the SVN repository of the open source 3D rendering engine, Ogre. The code is in the sample called DeferredShading.

Deferred Rendering In A Nutshell

Deferred rendering is an alternative approach to rendering 3d scenes. The classic rendering approach involves rendering each object and applying lighting passes to it. So, if an object is affected by 6 lights, it will be rendered 6 times, once for each light, in order to accumulate the effect of each light. This approach is often called "forward rendering".

Deferred rendering takes another approach: first of all of the objects render their "lighting related infomation" to a texture, called the G-Buffer. This includes their colours, normals, depths and any other info that might be relevant to calculating their final colour. Afterwards, the lights in the scene are rendered as geometry (sphere for point light, cone for spotlight and full screen quad for directional light), and they use the G-buffer to calculate the colour contribution of that light to that pixel.

The motive for using deferred rendering is mainly performance related – instead of having a worst case batch count of num_objects * num_lights (if all objects are affected by all lights), you have a fixed cost of num_objects + num_lights. There are other pros and cons of the system, but the purpose of this article is not to help decide whether deferred rendering should be used, but how to do it if selected.

The Problem

The main difficulty with implementing deferred rendering is that you have to do everything on your own. The regular rendering approach involves rendering each object directly to the output buffer (called 'forward rendering'). This means that all of the transform & lighting calculations for a single object happen in a single stage of the process. The graphics API that you are working with (DirectX, OpenGL etc) exposes many options for rendering objects with lights. This is often called the 'fixed function pipeline', where you have API calls that control the fashion in which an object is rendered. Since we are splitting up the rendering to two parts, we can not use these faculties at all, and have to re-implement the basic (and advanced) lighting models ourselves in shaders. Even shaders written for the forward pipeline won't be usable, since we use an intermediate layer (the G-Buffer). They will also need to be modified to write to the G-Buffer / read from the G-buffer (usually the first). In addition to that, the architecture of the rendering pipeline changes – objects are rendered regardless of lights, and then geometric representations of the light's affected area have to be rendered, lighting the scene. This is very different from the classic way – since when do we render lights?

The Goal

We would like to create a deferred rendering pipeline that is as unobtrusive as possible – we do not want the users of the engine to have to use it differently because of the way that its rendered, and we really don't want the artists to change the way they work just because we use a deferred renderer. So, we want an engine that can:
Interact with the game engine in the same way that the forward renderer does.
Use the same art and assets as the forward renderer and generate the same results
Be extended by users of the framework to still be as flexible as forward rendering.

Starting Point

This is not an article about starting a graphics engine from scratch. It assumes that you already have a somewhat high-level engine set up. The code that will be presented here is based on the open source Ogre 3D engine, but can be treated as semi-pseudo-code that can be re-implemented for other engines with similar concepts – materials, render order grouping and render-to-texture (including multiple render targets, referred to as MRT) pipeline control. Here is a short rundown of what would be required from a 3d rendering engine to build a deferred renderer on:

Material System: The system that stores all of the information that is required to render a single object type besides the geometry. Links to textures, alpha settings, shaders etc. are stored in an object's material. The common material hierarchy includes two levels:
Technique – When an object will be rendered, it will use exactly one of the techniques specified in the material. Multiple techniques exist to handle different hardware specs (If the hardware has shader support use technique A, if not fall back to technique B), different levels of detail (If object is close to camera use technique 'High', otherwise use 'Low'). In our case, we will create a new technique for objects that will get rendered into the G-buffer.
Pass – An actual render call. A technique is usually not more than a collection of passes. When an object is rendered with a technique, all of its passes are rendered. This is the scope at which the rendering related information is actually stored. Common objects have one pass, but more sophisticated objects (for example detail layers on the terrain or graffiti on top of the object) can have more.
Examples of material systems outside of Ogre are Nvidia's CGFX and Microsoft's HLSL FX.

Render Queues / Ordering System: When a scene full of objects is about to be rendered, who gets rendered when? That is the responsibility of this system. All engines need some control over render order since semi-transparent objects have to be rendered after the opaque ones in order to get the right output. Most engines will give you some control over this, as choosing the correct order can have visual and performance implications (less overdraw = less pixel shader stress = better performance, for example).

Full Scene / Post Processing Framework: This is probably the most sophisticated and least common of the three, but is still common. Some rendering effects, such as blur and ambient occlusion, require the entire scene to be rendered differently. We need the framework to support directives such as "Render a part of the scene to a texture", "Render a full screen quad", allowing us to control the rendering process from a high perspective.

In the Unreal Development Kit, this is called the "Unreal Post Process Manager". When working directly with OpenGL/ DirectX, you will have to write a layer like this on your own, or hardcode the pipeline (less recommended if creating an engine that is supposed to suit more than one game). Having a strong framework will also open up other possibilities like inferred lighting.

Generating the G-Buffer

So we know what we want to do, we can now start creating a deferred rendering framework on top of the engine. The problem of deferred rendering can be split up into two problems – creating the G-Buffer and lighting the scene using the G-Buffer. We will tackle both of them individually.

Deciding on a Texture Format

The first stage of the deferred rendering process is filling up a texture with intermediate data that allows us to light the scene later. So, the first question is, what data do we want? This is an important question – it is the anchor that ties both stages together, so they both have to synchronized with it. The choice has performance (memory requirements), visual quality (accuracy) and flexibility (what doesn't get into the G-Buffer is lost forever) implications.

We chose two FLOAT16_RGBA textures, essentially giving us eight 16 bit floating point data members. It's possible to use integer formats as well. The first one will contain the colour in RGB, specular intensity in A. The second one will contain the view-space-normal in RGB (we keep all 3 coordinates) and the (linear) depth in A.

Choosing a texture format is not a trivial decision, and has quite a big impact. The factors that come into play when choosing a texture format are:

GPU Memory consumption – This buffer will need to be created at the size of the viewport, for each viewport. For example, if the viewport is 1024x1024, you are paying one megabyte per byte of information per pixel. So, in our case, we have eight two-byte channels. That’s 16 bytes per pixel, so the G-buffer will require 16 megabytes of texture memory!
Flexibility – The G-Buffer is the ONLY link between the objects and the lighting process. Data that does not get into the G-buffer is lost forever. For example, we save enough information to light the scene using the standard model, but we do not know what object each pixel came from. Object-specific lighting options (for example, highlight the selected character) need an 'object ID' saved to the G-buffer as well.
Accuracy – More bits per data = more accuracy. Are 16 bits enough accuracy for depth? Maybe 8 bits per color channel would have been enough? This is usually a direct tradeoff with the memory factor.
Speed – Also a tradeoff with memory consumption. Some math tricks can be used to save memory. For example, the normal's Z coordinate can be recalculated from the X and Y coordinates if its sign is stored in a single bit (since X^2 + Y^2 + Z^2 = 1 => |Z| = Sqrt(1 – Y^2 – X^2) ), but those kind of calculations takes time.

In this case we chose the FLOAT16_RGBA pixel format mostly for simplicity. The information for the basic lighting model is there, and its easy to access. This can be changed later of course.

Preparing the objects for G-Buffer rendering

We assume that the material system has a scheming/profiling feature. This means that materials can specify multiple techniques that will be used in different scenarios (see 'Starting point' for futher explanation). We will use this system and define a new scheme – the G-Buffer scheme that will output the intermediate contents to the texture. Since we didn't want the artists to change the materials that they generate, the techniques for this scheme should be generated programatically (this can happen offline or at runtime) We will create this technique in all of the materials in our app. There are two ways to do this – either offline (load all the materials, add the GBuffer technique, save them back to disk) or online (when a material is loaded or about to be used, add the technique to its list). The lazy online approach was used in the Ogre sample, but all methods are possible.

How do we create the G-Buffer technique?

A) Inspect the classic technique

For each material (described in 'Starting Point'), we look at the technique that would have been used normally. For each pass in this technique, we check its properties and see what it does - does this pass have a texture? A normal map? Is it skinned? Transparent? Etc. The resulting classification should contain all the information required to build a GBuffer technique for an object. In some cases (like semi-transparent objects), it should also be able to flag that a certain object can not be deferred rendered, and it will be rendered normally later. In the Ogre sample, the function that does this is GBufferSchemeHandler::inspectPass, which returns a structure called PassProperties, which contains the fields needed to generate a g-buffer writing technique.

B) Generate the G-Buffer technique

After a pass has been inspected and understood, the next stage is to generate the G-Buffer-writing technique. Since the classification has all the information needed, it is possible to derive the material properties and create the matching shader to write to the MRT. A question that arises here is how should the shader be generated? There are two common approaches for this – the ubershader approach and the shader generation approach.

The first one involves writing a big shader beforehand with lots of preprocessor directives and compiling the right version of the shader during runtime, and the latter means generating shader code on the fly. For this task, I chose code generation as the option for two reasons:
There are quite a few options that affect G-Buffer rendering – diffuse textures, normal / specular / parallax maps, vertex colours, skinning etc. Writing a single ubershader for this will be very hard, because of the many preprocessor definitions. For example : how do you assign texture coordinate indices?
Debugging is comfortable because you see a simple shader that does exactly what a certain object needs. Here is an example of a shader pair generated (in this case, a normal map, a texture and a single predefined specularity level)
For example:

void ToGBufferVP( float4 iPosition : POSITION, float3 iNormal   : NORMAL, float2 iUV0 : TEXCOORD0, float3 iTangent : TANGENT0, out float4 oPosition : POSITION, out float3 oViewPos : TEXCOORD0, out float3 oNormal : TEXCOORD1, out float3 oTangent : TEXCOORD2, out float3 oBiNormal : TEXCOORD3, out float2 oUV0 : TEXCOORD4, uniform float4x4 cWorldViewProj, uniform float4x4 cWorldView ) { oPosition = mul(cWorldViewProj, iPosition); oNormal = mul(cWorldView, float4(iNormal,0)).xyz; oTangent = mul(cWorldView, float4(iTangent,0)).xyz; oBiNormal = cross(oNormal, oTangent); oViewPos = mul(cWorldView, iPosition).xyz; oUV0 = iUV0; } void ToGBufferFP(         float3 iViewPos : TEXCOORD0,         float3 iNormal   : TEXCOORD1,         float3 iTangent : TEXCOORD2,         float3 iBiNormal : TEXCOORD3,         float2 iUV0 : TEXCOORD4,         out float4 oColor0 : COLOR0,         out float4 oColor1 : COLOR1,         uniform sampler sNormalMap : register(s0),         uniform sampler sTex0 : register(s1),         uniform float4 cDiffuseColour,         uniform float cFarDistance,         uniform float cSpecularity         ) {         oColor0.rgb = tex2D(sTex0, iUV0);         oColor0.rgb *= cDiffuseColour.rgb;         oColor0.a = cSpecularity;         float3 texNormal = (tex2D(sNormalMap, iUV0)-0.5)*2;         float3x3 normalRotation = float3x3(iTangent, iBiNormal, iNormal);         oColor1.rgb = normalize(mul(texNormal, normalRotation));         oColor1.a = length(iViewPos) / cFarDistance; }

Remember to keep the shaders (uber or generated) synchronized with the G-Buffer format you decided on. Here is a screenshot from NVIDIA's PerfHUD tool showing the G-Buffer being built:

You can see the G-buffer textures (two of them) on the right side of the picture. Notice that nothing has been written to the back-buffer yet.

In the Ogre sample, the GBufferMaterialGenerator class takes care of this.

C) Postpone transparent / other none-deferred objects

The pass inspection tells us if a pass can be deferred or not. If not, we want to make sure we can render the none-deferred objects later. Again, material techniques come into play, and a NoGBuffer technique is introduced. When a none-deferred pass is detected, it is copied as-is to the NoGBuffer technique, allowing us to tell the framework to render all none-deferred objects once the deferred composition part is over. We will get to that later.

Allow overriding the automatic process

The shader generation makes it easier to manage assets in a deferred rendering environment since it generates the shaders from the fixed function parameters, but this is not always possible. You might have objects with specialized shaders that don't fit in any generic scheme, but still want to defer them as well. It is worthwhile to keep this option open, for better flexibility. The proposed framework does this already, as the programmatic material/shader generation only happens when an object without a technique for the GBuffer scheme defined. This means that if an object has a technique predefined, it will override the automatic process. This makes upkeep harder (you need to synchronize all manual shaders with the MRT format if you change it, for example) but is unavoidable in certain scenarios.

G-Buffer Generation Summary

Using these four stages, we added a hook in the system that when the GBuffer scheme is enabled, existing objects will be inspected, classified and have matching writing materials and shaders generated for them, allowing the existing pipeline to render itself to the G Buffer without asset modification. Some objects will be postponed and forward rendered later, some objects will be rendered using custom shaders – flexibility and compatibility remains, and the G-Buffer is created!

Lighting the Scene

We now have a prepared G-buffer with all the intermediate data we need to light the scene. Our next job is to render each light to the scene, calculating its contribution to the final image. This stage differs from the standard approach in a big way – we render lights! There is a question of who triggers the actual rendering. For that, we use the full scene / post processing framework (See description in 'Starting point'). In Ogre, it is called the compositor framework.

Compositing Scenes (a general idea)

The classic forward rendering pseudo code is:

for each visible object: for each light that affects object: render object with light contribution to main window

However, this is not always the case, even before deferred rendering. Many post-processing effects such as blurring the scene require rendering the scene to a texture, and then rendering the texture to the final output image using a pixel shader with a full screen quad. A simple motion blur pseudo code might be:

for each visible object: for each light that affects object: render object with light contribution to a texture blend texture with 'previous frame' texture to screen copy texture to 'previous frame' texture

A good scene composition framework will allow these kind of pipelines to be defined. Ogre does this with the compositor framework, which allows these kind of pipelines to be defined in scripts. For example, this is the GBuffer generating compositor:

compositor DeferredShading/GBuffer { technique { texture GBufferTex target_width target_height PF_FLOAT16_RGBA PF_FLOAT16_RGBA chain_scope target GBufferTex { input none pass clear { } shadows off material_scheme GBuffer pass render_scene { //These values are synchronized with the code first_render_queue 10 last_render_queue  79 } } } }

Even if you never heard of Ogre, this script should be self explanatory. The multi render texture is defined as two FLOAT16_RGBA textures, it is cleared and the scene is rendered to with the GBuffer material scheme every frame.

Back to our case…

This is exactly the time of custom control we need. We will now define a custom composition pass that will render light as geometry.

compositor DeferredShading/CompositeScene { technique { //Reference the main Gbuffer texture texture_ref GBuffer DeferredShading/GBuffer GBufferTex target_output { input none //We will dispatch the shadow texture rendering ourselves shadows off // render skies and other pre-gbuffer objects pass render_scene { first_render_queue 1 last_render_queue  9 } //Render the lights and their meshes pass render_custom DeferredLight { input 0 GBuffer 0 input 1 GBuffer 1 } pass render_scene { material_scheme NoGBuffer first_render_queue 10 last_render_queue 79 } } } }

Once a proper full scene composition framework is in place, this task becomes pretty simple to design. This compositor uses the result of the GBuffer compositor (see the texture_ref declaration), renders the skies (early render groups), then the light's geometry, then the objects that we skipped earlier (this is where the NoGBuffer scheme comes into place). The only thing left to do is to implement the deferred light composition pass. The main challenge is the fact that in deferred rendering, all the lighting calculations are your responsibility – including basic diffuse/specular lighting.

Once again, the question of shader management comes into play. This time around I actually chose to use the uber-shader approach, because the inputs don't change frequently, and there are less combinations to manage. Light type, and a shadow casting flag are enough for almost everything, so it doesn't create a mess. There are many papers on how the calculations are made, and the demo code is also a reference, although a basic one. The basic rule of thumb is that the G-buffer gives you view-space position and normal, which should be enough.

Creating the geometry representing the light is not a difficult task, creating a quad / sphere / cone programmatically is basic 3d geometry stuff. The attached demo contains code that does that. (See the GeomUtils class).

There are other minor issues – depending on your framework, you might need to reconstruct the depth buffer for the future objects to be able to depth test against the scene. Since the G-Buffer contains the depth as well, it is possible to rebuild the depth buffer from it, by rendering a quad with a pixel shader that reconstructs view space position, multiplies by the projection matrix and outputs the depth (remember – pixel shaders can output depth). In the framework, this is done by the ambient light, which is a 'fake light' that fills the scene with the ambient color and rebuilds depth, and is also rendered during the light composition pass. (Solutions to this problem are also already explained to good detail on the internet).

Here are some screenshots illustrating lighting of the scene using the G-Buffer :

As you can see, this time around the G-Buffer textures are used as input to the light geometry rendering, with the 3rd texture being the shadow map built for the light (it is built during the scene lighting stage to re-use the same shadow texture for all lights).

Here is a visualization of the scene after one light has rendered itself to the main buffer :

This light is a spotlight that casts shadows. Once all the lights render themselves in this fashion, the scene is fully lit!

Undeferred objects will be rendered normally afterwards – a render scene directive will be issued in the full scene rendering framework after the lighting takes place. This is a drawback since you will have to implement forward rendering techniques for them (which might lead to code duplication), but objects that can't be deferred are usually different from normal objects, so they would require special treatment anyways.

Additional post processing effects

Another advantage of deferred rendering, is that some advanced post-filters (like SSAO) require full-scene renders to get intermediate information about the scene and use it. If the G-Buffer contains this information, you can apply these effects without another full scene render; just make sure that your framework makes it easy to pass information (such as textures) from different render sequences to each other, and you are practically done. In the case of OGRE, the texture_ref directive (also used in the final scene composition) is all that is needed. The deferred rendering demo that accompanies this paper contains an SSAO postfilter, so you can switch it on and off and see the visual/performance impact.

Summary

My goal in this article was to get into some of the less intuitive details of deferred rendering implementations. I hope that this article will help anyone trying to implement a deferred rendering framework get to their target. Remember that there is full source code of a sample implementation as part of the Ogre SDK, so you can see it in action for full reference. For people with direct ogre experience, there is also an article that explains the Ogre usage in the demo. Good luck!

Cube Map Rendering Techniques for Direct3D10

Mon, 01 Feb 2010 09:32:33 +0000
I have recently been working on porting the Kourjet engine from Direct3D9 to Direct3D10. One of the things I wanted to improve was the rendering of point light shadow maps. In Direct3D9, you do not have much choice in rendering a cube map. The limitations of the API only allow rendering the six faces of the cube map in six distinct passes. However, new features in Direct3D10, such as geometry shaders and improvements to geometry instancing allow you to implement alternative algorithms for rendering cube maps in a single pass. This article will present three such algorithms and discuss their drawbacks and advantages.

Before we get to each algorithm, we will introduce some notations and conventions that will be used in the article. Our goal in the following algorithms is to render as efficiently as possible an arbitrary cube map located in an arbitrary location in the scene. The cube we are rendering to is aligned on the axes of the world basis (X, Y, Z).

It is assumed in this article that the reader is already familiar with the concept of cube mapping as well as the fundamentals of working with vertex and geometry shaders.

Algorithm 1 - Multi-Pass Rendering
Multi-Pass Rendering is the simplest algorithm and also the only one that can be used on Direct3D9-class hardware. The idea is to render the scene to each one of the cube faces separately. The (naive) algorithm could look like this:

BaseObjectSet = Select Objects to be Rendered on the Cube Map (1) For each face F of the cube map Set the Render Target corresponding to F Set the view and projection matrices corresponding to F Foreach Object O in BaseObjectSet     Render( O ) End End (1) is a step that allows you to select a subset of the scene according to the context of the cube map. For example, if you are rendering a cube map for shadow mapping, you will here select only the objects in the light range and casting shadows. If it is a reflection cube-map, you could include only objects close enough to the cube position.

If BaseObjectSet contains N objects, we will then render exactly 6N objects and make six render target switches.

To improve on this first version of the algorithm, we can notice that we could apply one more culling step: as we are rendering one face at a time, we can test the current object to see if it is in the current view frustum (constructed from the view/projection matrices for the current face). The improved algorithm would then look like this:

BaseObjectSet = Select Objects to be Rendered on the Cube Map (1) For each face F of the cube map Set the Render Target corresponding to F Set the view and projection matrices corresponding to F Compute the view frustum VF corresponding to F Foreach Object O in BaseObjectSet     If O is visible in VF Then         Render( O ) End End We can thus skip objects that are not visible on all the faces. An object of standard size such as a character or a car, will usually only be visible on two to three faces at one time. Thus, at worst we still render 6N objects but in practice we will render on average 2N to 4N objects (this of course depends on the exact set of objects in the scene). We however still have to switch the render target six times.

Algorithm 2 - Single-Pass Rendering Using Only a Geometry Shader
Direct3D10 comes with a sample program that renders a cube map in a single pass using a geometry shader. The idea is to pass six view/projection matrices to the geometry shader, and for each object primitive, generate six other primitives that are projected onto their respective cube face. We can thus take our previous algorithm and factorize some operations:

BaseObjectSet = Select Objects to be Rendered on the Cube Map (1) Set the six Render Targets corresponding to all the faces as an array Set the six view and projection matrices corresponding to all the faces as an array Foreach Object O in BaseObjectSet Render( O ) End As you can see, we got rid of the cube face loop. Before iterating the object set, we compute and set the needed variables for all faces at once (matrices and render targets). Similarly, there will be only one draw call per object. The number of API calls has been divided by six here compared to the worst-case scenario! Here is what the vertex shader and geometry shader could look like:

VS_OUTPUT VertexShader( VS_INPUT input ) { VS_OUTPUT output = (VS_OUTPUT) 0; output.position = mul( float4( input.position, 1), matWorld );  // Position in world space return output; } [maxvertexcount(18)] void GeometryShader( triangle VS_OUTPUT input[3], inout TriangleStream CubeMapStream ) { // For each face of the cube, create a new triangle [unroll] for( int f = 0; f < 6; ++f ) {     GS_OUTPUT output = (GS_OUTPUT) 0;                      // Assign triangle to the render target corresponding to this cube face     output.RTIndex = f;                      // For each vertex of the triangle, compute screen space position and pass distance     [unroll]     for( int v = 0; v < 3; ++v )     {         output.position = mul( input[v].position, matViewProj[f] );         CubeMapStream.Append( output );     }                      // New triangle     CubeMapStream.RestartStrip(); } } The vertex shader simply transforms the vertices from object space to world space. Most of the work is done in the geometry shader. In fact our face loop from the first algorithm can be found here: we have traded CPU cycles for GPU cycles. One problem that seems obvious in this naive implementation is that there is no culling on a face basis. Which means that the geometry shader will always be outputting six times more vertices.

One way to go around this problem is to test the object against each of the faces' respective frusta and set a flag to be used in the geometry shader to skip generation of some primitives. We can then change the algorithm as follow:

BaseObjectSet = Select Objects to be Rendered on the Cube Map (1) Set the six Render Targets corresponding to all the faces as an array Set the six view and projection matrices corresponding to all the faces as an array Compute the six view frusta corresponding to all the faces Foreach Object O in BaseObjectSet For each face F of the cube map     Test O with the F's frustum     If O is visible, set the geometry shader flag for F (cubeFaceWriteFlag |= 1 << F) End Render( O ) End [maxvertexcount(18)] void GeometryShader( triangle VS_OUTPUT input[3], inout TriangleStream CubeMapStream ) { // For each face of the cube, create a new triangle [unroll] for( int f = 0; f < 6; ++f ) {     // Only output primitives to the required cube faces     [branch]     if ( ( cubeFaceWriteFlag & ( 1 << f ) ) != 0 )     {         GS_OUTPUT output = (GS_OUTPUT) 0;                          // Assign triangle to the RT corresponding to this cube face         output.RTIndex = f;                          // For each vertex of the triangle, compute screen space position and pass distance         [unroll]         for( int v = 0; v < 3; ++v )         {             output.position = mul( input[v].position, matViewProj[f] );             CubeMapStream.Append( output );         }                          // New triangle         CubeMapStream.RestartStrip();     } } } Although we have managed to skip generation of some primitives, the geometry shader is still not optimal. Indeed, the geometry shaders are usually optimized for generation of at most four vertices whereas we are generating at most 18.

Algorithm 3 - Single-Pass Rendering Using a Geometry Shader and Geometry Instancing
Geometry instancing is a concept that has already been around in Direct3D9. However, Direct3D10 and its geometry shaders open some new ways of using it. In Direct3D9, one would pass a vertex buffer and a buffer containing instance data (for example six transformation matrices) and the API will, in a single draw call, render the vertex buffer six times, each with different instance data. Very nice, but useless for us when we need to render to a cube map: there is no way to select dynamically the render target that will be used for a given instance.

This is where the geometry shaders come into picture. It is indeed possible as seen in the previous algorithm to indicate per primitive to which render target it should be rendered. The idea will then be to have an instance buffer dynamically filled with the face indices on which the object is visible. So here we are with a first version of our third algorithm:

Select Objects to be Rendered on the Cube Map (1) Set the six Render Targets corresponding to all the faces as an array Set the six view and projection matrices corresponding to all the faces as an array Compute the six view frusta corresponding to all the faces Foreach Object O in BaseObjectSet Reset Dynamic Instance Buffer For each face F of the cube map     Test O with the F's frustum     If O is visible, append F to the instance buffer End Render six instances of O using our dynamic instance buffer End The algorithm is pretty similar to the previous one: we manage to minimise the number of API calls by rendering to all the cube faces at once. By testing against the face frusta, we only append to the instance buffer the relevant faces, thus not instancing any invisible geometry (as in the previous algorithm). Even though this looks pretty smooth in our pseudo algorithm, there are some hidden complications.

In Direct3D, as we are free to organise the vertex/instance data as we wish, we must tell the API (using vertex layouts) how data we are providing is organised. Usually the engine will store this per mesh and switch the vertex layout when drawing the object. Indeed, we want to keep the freedom of being able to specify any vertex layout (maybe one for characters, one for terrain, one for sky, etc.). In this algorithm, we are appending instance data and thus have to append some information to the existing vertex layout. In the Kourjet Engine for example, the new vertex layouts are computed on the fly the first time we draw an object and kept in a cache belonging to the render path.

The vertex shader and geometry shader are pretty simple once we got the vertex layout right. We basically use the vertex shader to do all transformations and the geometry shader simply to direct the primitive to the appropriate render target. Here is the source code:

struct VS_INPUT { float3 position : POSITION; uint   cubeFace : CUBEFACE; }; struct VS_OUTPUT { float4 position : SV_POSITION; // Position in screen space uint   cubeFace : CUBEFACE; }; struct GS_OUTPUT { float4 position : SV_POSITION; // Position in screen space uint   RTIndex  : SV_RenderTargetArrayIndex; }; VS_OUTPUT VertexShader( VS_INPUT input ) { VS_OUTPUT output = (VS_OUTPUT) 0; float4 P = mul( float4( input.position, 1), matWorld ); // Position in world space output.position = mul( P, matViewProj[ input.cubeFace ] ); output.cubeFace = input.cubeFace; return output; } [maxvertexcount(3)] void GeometryShader( triangle VS_OUTPUT input[3], inout TriangleStream CubeMapStream ) { GS_OUTPUT output = (GS_OUTPUT) 0; // Assign triangle to the render target corresponding to this cube face output.RTIndex = input[0].cubeFace; // For each vertex of the triangle [unroll] for( int v = 0; v < 3; ++v ) {     output.position = input[v].position;     CubeMapStream.Append( output ); } // New triangle CubeMapStream.RestartStrip(); }
Conclusion

The algorithm pros and cons
    Pros Cons   Algorithm 1
Works on Direct3D9 hardware
Can do visibility culling on each face's frustum

Need to switch render targets six times
At worst six draw calls per object
    Algorithm 2
One render target switch
One draw call per object

Geometry shader could be less than optimal depending on the hardware and the vertex format ^[1]
No real culling per face (vertex shader is still run six times)
Direct3D10 only
    Algorithm 3
One render target switch
One draw call per object
Can do visibility culling on each face's frustum (instancing only geometry if visible on a face)

Must maintain a cache of duplicated vertex layouts for instancing
Direct3D10 only

Source code and implementation
This article does not come with any particular source code attached. However implementation of these algorithms applied to rendering of shadow maps can be found in the Kourjet Engine source code, and in particular the following files can be of interest:

Algorithm 1 and its corresponding effect file
Algorithm 2 and its corresponding effect file
Algorithm 3 and its corresponding effect file
Reference articles & other links

The blog posts about the implementation of these algorithms in Kourjet for shadow mapping (part 1, part 2 and part 3)
In this paper on GPU programming, NVidia gives some recommandations for Geometry Shaders, and in particular in section 4.6, page 40, they say:

The performance of a GS is inversely proportional to the output size (in scalars) declared in the Geometry Shader, which is the product of the vertex size and the number of vertices (maxvertexcount). This performance degradation however occurs at particular output sizes, and is not smooth. For example, on a GeForce 8800 GTX a GS that outputs at most 1 to 20 scalars in total would run at peak performance, whereas a GS would run at 50% of peak performance if the total maximum output was declared to be between 27 and 40 scalars.

More concretely, if your vertex declaration is:

float3 pos: POSITION; float2 tex: TEXCOORD;
Each vertex is 5 scalar attributes in size.

If the GS defines a maxvertexcount of 4 vertices, then it will run at full speed (20 scalar attributes).
But if you increase the vertex size by adding a float3 normal, the number of scalars will increase by 4 * 3, putting your total at 32 scalar attributes. This will result in the GS running at 50% peak performance.

Walls and Shadows in 2D

Tue, 10 Nov 2009 04:44:27 +0000
This article will explore a technique for efficiently computing visible and lit regions in a 2D scene.

The algorithms presented are intended for use with large maps, and where computation time of some appreciable fraction of a second is tolerable. There are faster algorithms, but they involve solving the intersection and union of arbitrary polygons; the reasons for rejecting this approach are discussed below.

Consider a region, such as a cave, dungeon or city, viewed from above, such that floors are shown by areas and walls by lines. Floors are always enclosed by some set of walls. Put another way, all maps are bounded by a continuous wall. We will ignore other map features; for our purposes, the only “important” features are walls and floors.

Lights exist at points, and have a given area of effect, defined by a radius. Everything inside this radius that faces the light is considered lit. Lights can have overlapping areas. Walls cast shadows, making the walls and floors behind them unlit (there is no provision for light attenuating – getting dimmer over distance – things are either lit or unlit).

The observer also exists at a point. The problem is to quickly discover what he can see, where the area seen is defined as any point he has line of sight to (that is, those walls and floors not occluded by other walls), that are also lit. For efficiency reasons, the observer’s vision also has a maximum radius. Anything outside that radius is not visible, even if it is lit and unobstructed. The algorithm can function without this limit, the limitation is an optimization.

In addition to knowing what areas are currently visible, it is desirable to know what areas have been previously visible, and to be able to display those areas differently.

An example is given in Figure 1. The observer is the red hollow circle, and standing at the intersection of three passageways. White solid circles are light (the observer, in this case, is carrying one).

Dark green is used for walls and floors that have not been seen yet. Blue walls and medium grey floors mark areas that have been seen previously. White walls and light grey floor denote what’s currently visible. Note that while the observer’s own light doesn’t reach all the way down that west passageway, there’s a light to the south that illuminates part of it, giving a disjoint area of visibility. To the southeast, there’s another light that helps light that southeast corridor, so the whole corridor is visible. Finally, to the northeast, there’s a small window in the wall, allowing the observer to illuminate, and see, a small amount of the room to the east.

Figure 1 - A simple example of a map In this technique, the floor is broken into small fragments (triangles in this case, because they are easy to render), which serve as the algorithm’s basic unit of area. When determining what part of the floor is visible, I’m really asking what set of triangles is visible. If the centroid (average of the vertices) of a triangle is visible, the whole triangle is considered visible. This has the side effect of making the bounds of the visibility area slightly jagged, as you can see in the disjoint area down the west corridor. In my own application, this is acceptable, but blending of adjacent triangles could be done in order to get a smooth gradient between visible and invisible areas.

Because walls divide up floor area and walls can run at all sorts of angles, cutting the floor into small triangles often results in additional, smaller, triangles being generated. All told, a map can contain millions of floor triangles and hundreds of walls. For the curious, figure 2 shows the triangles generated for a small section of the map:

Figure 2 - You are in a maze of twisty triangles, all largely the same For lighting (and seeing) the walls themselves, I do something similar – walls are broken into segments, and the center point of each segment is checked to see if it is visible. If it is visible, that whole segment is visible. For this reason, segments are kept short and walls can generate thousands of them in a large map.

Because of this, when it comes time to determine which parts of walls and floors are not visible it may be necessary to evaluate millions of points for the floor and thousands of points for wall segments. Conceptually, they all need to be evaluated against every wall to determine if line of sight exists from the observer, and that process has to be repeated for each light as well.

Clearly, a brute force approach will not work in reasonable time. The goal is to move the observer to a new point, or move a light to a new point (often both, since the observer often carries a light), and know as quickly as possible what areas of floor and segments of walls can be seen. Comparing possibly millions of points against hundreds or thousands of walls and doing a line of sight calculation – essentially calculating the intersection of two line segments, for one for a line of vision and one for a wall – isn’t acceptably fast.

It turns out that lighting and vision can be handled by the same algorithm, since they are both occluded by walls in the same way. They can both be represented by casting rays out from a given point, and stopping the rays when they hit a wall. If there’s no wall along that ray, the ray is cut short by a distance limit instead (this illustrates another difficulty with using polygon intersections – polygons don’t have curved sides, and approximating them with short straight lines increases the cost of the intersection test).

Figure 3 - A "polygon" of light, with curved parts Since there can be multiple lights, unions of polygons would be required:

Figure 4 - Union of two lit areas Since we’ve defined visibility as areas that are both lit and within line of sight of the observer, the intersection of polygons representing lit areas (itself a union) and the polygon representing the area of sight represent the visible areas. Figure 5 repeats the original example, with yellow lines roughly delineating the lit area, and red lines bounding the area of sight. Areas within both are the visible areas.

Figure 5 - Vision and Light compute Visibility The result of the intersection is a (potentially empty) set of polygons. To maintain a history of what’s been visible, a union of the previously visible areas and the currently visible area is also required. The combination can create polygons with holes, and for complex maps, a large number of sides. Figure 6 shows the result of a wall, a number of square columns, and a short walk along the north side of the wall by the observer, carrying a light. The union of previously visible areas is shown in darker grey¹.

Figure 6 - Columns, a Wall, and A Messy Polygon Union Doing polygon union and intersection is complex. Naïve implementations of these algorithms run into problems with boundary conditions and, in complex cases, floating point accuracy. There are packages available that solve these problems, and deal elegantly with disjoint polygons and holes, but they are available under restricted license. I wanted an unencumbered solution, and was willing to trade off some amount of runtime to get it.

But a brute force computation of every floor point vs. every obstructing wall is unacceptable. What is needed is an efficient way to evaluate the many floor and wall points for visibility.

The basic approach amounts to describing the shadows cast by walls. Since each point in the floor has to be tested against these shadows, the algorithms focus on making it as inexpensive as possible to determine if any given point is within a shadow, without loss of accuracy. Since “as inexpensive as possible” is still too expensive, given the sheer number of points to consider, the approach also includes determining which walls are already covered by other walls (in effect, we work to discard walls which don’t change anything we care about.) Where walls cannot be discarded, the algorithm attempts to determine what parts of walls contribute to meaningful shadows, and which parts are irrelevant. Note that while I talk about shadows here, everything also applies to occluding the observer’s view; since, as noted, a wall stops a line of sight in exactly the same way that it cuts short a ray of light. Finally, I discuss the critical optimizations that make the approach fast enough to use on large and complex maps.

We start by making all points are relative to the location of the light in question; in other words everything is translated so that the light is at the origin. This translation drops terms out of many formulas and provides significant savings.

The first step is to identify when a wall casts a shadow over a point – efficiently. The simplest way to do this is to take the endpoints of a wall and arrange them so that the first point is clockwise from the second point, as seen from the origin. If they are counterclockwise, swap the points. If they are collinear with the light, throw out that wall, because it can’t cast a shadow. I will refer to the endpoints as the left and right endpoints, with the understanding that this is from the perspective of the origin.

A 2D cross product³ (from the origin, to the start point, to the end point), reveals both the edge-on case and the clockwise or counterclockwise winding of the end points.

Here’s an example:

Figure 7 - Walls that do, and don't, matter In figure 7, line C is edge on to the light and casts no shadow of its own, so it gets dropped. B isn’t edge on, so we keep it, with BA as its right endpoint and BC as the left (as seen from the light at the origin). D is also a wall that matters, and DC becomes the right endpoint, while DE becomes the left. In some applications, in which walls always form continuous loops, A and E can also be dropped because they represent hidden surfaces. The same rules that apply to 3D surface removal apply here – in a closed figure, walls that face away from the origin can be dropped without harm. However, this algorithm works even if walls don’t form closed figures.

Now, cast a ray out from the origin though the right endpoint of a given wall – use B as an example, and cast a ray through BA. Note that any point that is in shadow happens to be to the left of this line (so are many other points, but the point is that all the shadowed ones are). Calculating “to the left of” is cheap: it’s the 2D cross product from the origin, to the right endpoint of the wall, to the point in question; for example, point X in Figure 7. The cross product gives a value that’s negative on one side of the wall, positive on the other side, and zero if the point in question is on the line.

Repeat for the ray from the origin to the left end of the line segment, BC. All shadowed points are to the right of this line, which is determined by another 2D cross product. All that remains is to determine if the point is behind the wall or in front of the wall, with respect to the center. This is yet another 2D cross product, this time from the wall’s right endpoint, to the left end point, to the point in question. Point X in Figure 7 would pass all three of these tests, so it is in B’s shadow.

All told, at most three cross products (six multiplies and three subtracts, and three comparisons with 0), tell if a point is shadowed by a wall. In many cases, a single cross product will prove that a point is not shadowed by a wall. But that still leaves the problem of comparing many, many thousands of points against hundreds of walls.

Having established an algorithm to test a point against a wall, we now need to find ways to minimize how often we have to use it. Any wall we can cull results in hundreds of thousands of fewer operations! So a first pass at culling is simply to remove any wall which is outside the radius of the light by creating a bounding square around the origin with the “radius” of the light, and a bounding rectangle from each wall’s endpoints. If these don’t overlap, that wall can’t affect lighting, and is discarded.

Figure 8 - Using rectangles and overlap to discard walls In Figure 8, F’s rectangle doesn’t overlap the light’s rectangle, so F gets discarded. H and G overlap and so are kept – H is a mistake because it’s not really in the circle of light that matters, but this is a very cheap test that discards most walls in a large map very quickly, and that’s what we want for now. In applications where walls form loops or tight groups, the entire set can be given a single bounding rectangle, allowing whole groups of walls to be culled by a single overlap test.

Whatever walls are left might cast shadows. For each, we calculate the squared distance between the origin and the nearest part of the wall. This is slightly messy, since “nearest” could be either endpoint, or a point somewhere between. Given these distances, sort the list of walls so that the closest rise to the top. In the case of ties, there is generally some advantage in letting the longer wall sort closer to the top. This will usually put the walls casting the largest shadows near the top of the list. This helps performance considerably, but nothing breaks if the list isn’t perfectly sorted. In figure 8, G would be judged closer, by virtue of the northernmost endpoint. H’s closest point, near the middle of H, is further off.

Once we’ve dropped all the obviously uninvolved walls and sorted the rest by distance, it’s time to walk through the list of walls, adding them to the (initially empty) set of walls that cast shadows. If a wall turns out to be occluded by other walls in this phase, we cull it. Usually, anyway - in the interest of speed, the algorithm settles for discarding most such walls, but can miss cases. In practice, it misses few cases, so I have not troubled to improve this phase’s culling method.

To explain how this culling is done, we must introduce some concepts. Each wall generates a pair of rays, both starting from the light (origin) and one through each endpoint. As noted before, points to the right of the “left” ray, and to the left of the “right” ray, bound the possible shadow. However, some of that area might be already shadowed by another wall – one wall can partially cover another. In fact, all of that area might be shadowed by other walls – the current wall might be totally occluded. If it’s only partially occluded, what we want to do is “narrow” its rays, pushing the left ray further right and/or the right ray further left, to account for the existing, known shadows in this area. The reason for this is that we don’t want to compare points against multiple walls if we don’t have to, and we have to compare a point against any wall if it lies between the wall’s left and right rays. The narrower that angle becomes, the fewer points have to be checked against that wall.

So when we take in a new wall, the first thing we do is look at the two line segments between the origin and the wall’s endpoints (each in turn). If that line segment intersects a wall we already accepted, then the intersected wall casts a shadow we probably care about in regard to the current wall. The interesting point here is that we don’t care where the intersection occurs. An example will show why:

Figure 9 - Intersecting origin-to-endpoint with other walls Assume the current wall, W, is half-hidden by an already accepted, closer one, S. Assume that it’s the left half of W that’s covered by the nearer wall, as in the example above. The way we discover this is by checking the line segment from the origin to left endpoint of W, against the already-accepted walls. It can intersect several. Once we find an intersection, we know immediately that the intersected wall, S, is going to cover at least part of W, and on the left side (ignore the case where the S’s left endpoint and S’s right endpoint are collinear with the origin – it doesn’t change anything). Notice that we don’t care where S gets intersected.

So what do we do? We replace W’s left endpoint ray with S’s right endpoint ray. In effect, we push the left endpoint ray of W to the right. Having done that, we check to see if we’ve pushed it so far to the right that it is now at or past W’s own right ray. If so, S completely occludes W and we discard W immediately. If not, we’ve made W “narrower”. In our example, W survived and its shadow got narrower, as marked in grey.

Figure 10 - Pushing W's left ray We keep doing this, looking for other walls to play the part of S that intersect W’s left or right end rays. When they do, we update (“nudge”) W’s left (or right) ray by replacing it with S’s right (or left) ray. If the endpoint rays of W meet or cross during this “nudging” process, W is discarded.

Figure 11 - Losing W In the example above, S has pushed W’s left ray, and J has pushed W’s right ray. A 2D cross product tells us that left and right rays have gotten past each other in the process, so W is judged to be completely occluded, and gets dropped. Otherwise, it survives, with a potentially much narrowed shadow, and it adds it to the set of kept walls. Note that if J had been longer to the left, it could have intersected both W’s red lines, and it would have pushed both W’s left and right rays by itself, forcing them to cross. This makes sense; it would have completely occluded W all by itself and we’d expect it to cause W to drop out.

Note that this algorithm doesn’t notice the case where a short, close wall casts a shadow over the middle of a long wall a little further off. In this case, both walls end up passing this check, and the long wall doesn’t get its rays changed (the only way to do that would be to split the long wall into two pieces and narrow them individually). This isn’t much of a problem in practice, because when it comes time to check points, we will again check the walls in order of distance from the origin, so the smaller, closer wall is likely to be checked first. Points it shadows won’t have to be checked again, so for those points the longer wall never needs to be checked at all. There are unusual cases where points do end up in redundant checks, but they are unusual enough not to be much of a runtime problem.

As we work through the list of candidates, we are generally working outward from the origin, so it’s not uncommon for more and more walls to end up discarded because they are completely occluded. This helps keep this part of the algorithm quick.

It remains to find a good way to detect intersections of line segments. We don’t want round-off problems (this might report false intersections or miss real ones, causing annoying issues), and we don’t care where the intersection itself actually occurs. It turns out that a reasonable way to do this is to take W’s two points, and each potential S’s two points, and arrange them into a quadrilateral. We take the 4 points in this order: S’s left, W’s left, S’s right, W’s right. If the line segments cross, the four points in order form a convex polygon. If they don’t, it is concave. An example serves:

Figure 12 - Detecting intersections R and S cross, so the (green) polygon formed by the 4 endpoints, is a convex kite shape. R and L don’t cross, so the resulting (orange) polygon isn’t convex; it’s not even simple.

It turns out that there is a fair amount of misinformation about testing for convex polygons on the ‘Net. Talk about counting changes in sign sounds interesting (and cheap), but I’ve yet to see an implementation of this that works in all cases, including horizontal and vertical lines. What I’ve ended up with is more expensive but gets all cases, even if two points in the quad happen to be coincident. I calculate the 4 2D cross products, going around the quad (in either order). If they are all negative OR they are all positive, it’s convex. Anything else is concave or worse. While not cheap (up to 8 multiplies and quite a few subtracts), we can stop as soon as we get a difference in sign. On hardware that can do floating subtracts in parallel, this is not too bad in cost.

By itself, that’s enough to discard unneeded walls in most cases, and minimize the scope of influence of the surviving walls. Just with what we have, it’s possible to determine if points are occluded by walls. But we’d still like it faster. We always want things faster, that’s why we buy computers.

Making it faster
Make sure all we just discussed makes sense to you, because we’re about to add some complications. There are four optimizations that can be applied to all this, unrelated to each other. 1. In my maps, walls (except doors) always touch one other wall at their endpoint. (They are really wall surfaces - just as in a 3D model, all the surface polygons are just that, surfaces, always touching neighbor surfaces along edges.) This leads to an optimization, though it is a little fussy to apply. An example serves best.

Imagine you’re a light at the origin, and over at x=5 there’s a wall stretching from (5,0) to (5,5), running north and south. It casts an obvious shadow to the east. We’ll call that wall A. But imagine that at (5,5) there’s another wall, running to (4,6), diagonally to the northwest. Call it B.

Figure 13 - Extending A's influence A has the usual left and right rays: the right ray passes through (0,5), the left through (5.5). B has its own rays, with a right ray at (5,5) and a left ray at (4.6).

Between them, they cast a single, joined shadow, wider than the shadows they cast individually. The shape of the shadow is complicated, but it’s worth noticing that for any point behind the line of A (that is, with x > 5, noted in green), that point is in shadow if it is between A’s right ray and B’s left ray. This is because B extends A, and (important point) it extends it by turning somewhat towards the light, not away. (It would also work if A and B were collinear, but in my system, collinear walls that share an endpoint become a single wall). There are also points B shadows that have X < 5, but when we are just considering A, it’s fair to say that we can apply B’s left bound instead of A’s left bound when asking about points that are behind the line A makes. A’s ability to screen things, given in light grey, has effectively been extended by the dark grey area.

I don’t take advantage of this when it comes to considering which points are in shadow, because all it does is increase the number of points that are candidates for testing for any given wall, and that doesn’t help. However, I do take advantage of this when determining what walls occlude other walls.

I do this by keeping two sets of vector pairs for each wall. The ones I’ve been calling left and right are the “inner” vectors, named because they tend to move inward, and their goal is to get pushed closer together by other walls, ideally to cross. But there is also a pair of left and right vectors I call the outer pair. They start at the endpoints of the wall like the inner ones do, but they try to grow outward. They grow outward when 1) I find the wall that shares an endpoint and 2) this other wall (in my scheme there can only be one) does not bend away from the light. This is an easy check to make – it’s another 2D cross product, from A’s left endpoint, to A’s right, to B’s right. If that comes up counterclockwise, A’s outer right vector gets replaced by a copy of B’s outer right vector (as long as that’s an improvement – it’s important to check that you’re really pushing the outer right vector more to the right.)

And note the trick affects both A and B. If B extends A’s right outer vector, then A is a great candidate for extending B’s left outer vector.

Applied carefully, the extra reach this gives walls helps discard distant walls very quickly in cases where there are long runs of joined walls. I find that for most maps, the difference this makes is not large, and given the work I put into getting it right, I might not have done it if I’d realized how little it helps most maps.

2. When considering points, it’s important not to waste time testing any given point against a wall that can’t possibly shadow it. Each point, after all, has three tests it has to pass, per wall. Is it to the right of the left vector, to the left of the right vector, and is it behind the line of the wall. On average, half of all points are going to pass that first test, for most walls. That means that a fair amount of the time, the second test is going to be needed for points that are not remotely candidates. And given huge numbers of points, that’s unacceptable.

Figure 14 - The futility of any one test Here, P is to the right of the left endpoint-origin line, so it’s a candidate for being shadowed by the wall. But so are Q and R, and they clearly aren’t going to be shadowed. It would be helpful, then, if we could only test a point against the walls that have a good shot of shadowing it.

Trigonometric solutions suggest themselves, but trig functions are much too expensive.

What I do is create a list of “aggregate wedges.” Each wall’s shadow, called a wedge, is compared with the other wedges of the other walls. If they overlap, they are added to the same set of wedges, and I keep track of the leftmost and rightmost ray among everything in the same set.

Figure 15 - Creating sets of shaodw wedges Of course, if you’re in the square room without a door (never mind how you got in), you end up with all the walls in the same set, and the rays that bound the set end up “enclosing” every point on the map! So this trick is useless in these kinds of maps. But in maps of towns, with many freestanding buildings and hence many independent wedges, you can often get whole groups of walls into a number of disjoint sets, and since each set has an enclosing pair of rays that covers all the walls in the set, you can test any given point against the set’s rays: if it’s not between them, you don’t have to test any of the walls in that set.

This sounds pretty, but it can be maddening to get right. You can have two independent sets, and then find a wall that belongs in them both, effectively merging two sets into one big one.

Figure 16 - Combining sets You end up doing a certain amount of set management. I have a sloppy and dirty way to do this which is reasonably fast, but it’s not pretty.

Another difficulty is knowing when a wall’s vectors overlap an existing set’s bounds. There are several cases. A wall’s vectors could be completely inside the set’s bounds, in which case the wall just becomes a member of the set and the set’s bounds don’t change. Or it can overlap to the left, joining the set and pushing the set’s left vector. Or it can overlap on the right. Or it can overlap on both sides, pushing both set vectors outward. Or it can be disjoint to that set. Keep in mind that a set can have rays at an obtuse angle, enclosing quite a bit of area. It’s surprisingly hard to diagnose all these cases properly. The algorithm is difficult to describe, but the source code is at the end of this document.

When all is said and done, this optimization is very worthwhile for some classes of maps that would otherwise hammer the point test algorithm with a lot of pointless tests. But it was not much fun to write.

3. This is my favorite optimization because it’s simple and it tends to help a great deal in certain, specific cases. When I’m gathering walls, I’m computing distance (actually, distance squared – no square roots were used in any of this) between the light and each wall. Along the way I keep track of the minimum such distance I see. If, for example, the closest distance to any wall is 50 units, then any point that is closer to the origin than 50 units can never be in shadow and doesn’t have to be checked against any wall. (Of course, if the light’s actual radius of effect is smaller than this distance, that value is used instead).

Figure 17 - Identifying points unaffected by walls If the light is close to a wall, this optimization saves very little. If it’s not, this can put hundreds or even thousands of points into “lit” status very quickly indeed. In order to avoid edge cases, I subtract a small amount from the minimum wall distance after I compute it, so there’s no question of points right at a wall being considered lit unduly.

4. This is my second favorite optimization because it’s simple, dirty and very effective. When I generate floor triangles, I use an algorithm that more or less sweeps though in strips. Because triangles are small, and generated such some neighboring triangles are adjacent in the list of triangles, odds are often very high that if a triangle is shadowed by a wall, the next triangle to consider is going to be shadowed by the same wall.

Figure 18 - Neighbors often suffer the same fate So the best optimization of all is to remember which wall last shadowed a triangle, and test the next triangle against that one first, always. After all, if a triangle is shadowed by any wall, it doesn’t have to be tested against any other; we don’t care about finding the “best” shadow, we just want to quickly find one that works. If the light happens to be close to a wall (which ruins optimization 3), this one can be very powerful. W, here, is likely to shadow about half the map.

One final trick – not exactly an optimization – has to do with the fact that this code runs on a dual core processor. I cut the triangle list in about half and give each half to a separate thread to run though. (Each thread has its own copy of the pointer used in optimization 4, so they don’t interfere with each other in any way.) This trick doesn’t always cut the runtime in half – it’s not uncommon for one thread to get saddled with most or all of the cases that require more checks – but it helps. Other speedups involve not using STL or boost, and sticking to old fashion arrays of data structures – heresy in some C++ circles, but the speed gains are worth any purist’s complaints.

What’s left is trivial. Each floor triangle, and short wall segment, have a set of bits, one for each possible light, and one for the viewer. If the object is not in shadow, the appropriate light’s bit is set in that object. If any such bit is set, the object is lit. There is also a bit for the observer, which as noted uses the exact same algorithm. If the object is lit, that algorithm is then run for that point for the observer, and if it comes up un-occluded, it is marked visible (and also marked “was once visible”, because I need to keep a history of what has already been seen). Moving a light is a matter of clearing all bits that correspond to that light in all the objects, and then “recasting” the light from its new location. In my application, most lights don’t move frequently, so not much of that happens.

My favorite acid test at the moment is an open town map with 3 million floor triangles and almost 7000 walls. With ambient light turned on (which means everything is automatically considered lit, so everything has to be checked for “visible to observer”), and a vision limit of 400 units (so about 336,000 possible triangles in range), my worst case compute times are about 0.75 seconds on a dual core 2Ghz Intel processor. Typical times are a more acceptable 0.4 sec or so. Kinder maps (underground caves, tight packed cities, open fields with few obstructions) manage times of well under 0.1 sec.

Some examples of my implementation:

Figure 19 - Room in underground city
Figure 20 - Ambient light, buildings and rock outcrops
Figure 21 - Limited lights, windows and doors Notice that lights in buildings shine out windows, and also enable peeks inside the building from outside, forming a number of disjoint areas of visibility.

  ¹To give a sense of the algorithm’s performance, computing the currently visible area in figure 5, and combining it with the previous area, took 0.003 seconds on a dual core 2Ghz processor. However, keep mind that figure 5 represents very small and simple map (containing less than 100,000 triangles).
²GPC, PolyBoolean and LEDA are among these packages. Some discussion of the problems of runtime, runspace, and accuracy can be found at http://www.complex-a5.ru/polyboolean/comp.html. Boost’s GTL is promising, but it forces the user of integers for coordinates, and the implementation is still evolving.
³Formally speaking, there isn’t a cross product defined in 2D; they only work in 3 or 7 dimensions. What’s referred to here is the z value computed as part of a 3D cross product, according to u.x * v.y - v.x * u.y.
Code Listing
What follows gives the general sense of the algorithms’ implementation. Do note that the code is not compiler-ready: support classes like Point, Wall, WallSeg and SimpleTriangle are not provided, but their implementation is reasonably obvious.

This code is released freely and for any purpose, commercial or private – it’s free and I don’t care what happens, nor do I need to be told. It is also without any warranty or promise of fitness, obviously. It works in my application as far as I know, and with some adjustment, may work in yours. The comments will show some of the battles that occurred in getting it to work. The code may contain optimizations I didn’t discuss above.

  enum Clockness {Straight, Clockwise, Counterclockwise}; enum Facing {Colinear, Inside, Outside}; static inline Clockness clocknessOrigin(const Point& p1, const Point& P2) { const float a = p1.x * P2.y - P2.x * p1.y; if (a > 0)        return Counterclockwise; // aka left if (a < 0)        return Clockwise; return Straight; } static inline bool clocknessOriginIsCounterClockwise(const Point& p1, const Point& P2) { return p1.x * P2.y - P2.x * p1.y > 0; } static inline bool clocknessOriginIsClockwise(const Point& p1, const Point& P2) { return p1.x * P2.y - P2.x * p1.y < 0; } class LineSegment { public: Point begin_; Point vector_; //begin_+vector_ is the end point inline LineSegment(const Point& begin, const Point& end)     : begin_(begin), vector_(end - begin) {} inline const Point& begin() const {return begin_;} inline Point end() const {return begin_ + vector_;} inline LineSegment(){} //We don't care *where* they intersect and we want to avoid divides and round off surprises. //So we don't attempt to solve the equations and check bounds. //We form a quadilateral with AB and CD, in ACBD order. This is a convex kite shape if the // segments cross. Anything else isn't a convex shape. If endpoints touch, we get a triangle, //which will be declared convex, which works for us. //Tripe about changes in sign in deltas at vertex didn't work. //life improves if a faster way is found to do this, but it has to be accurate. bool doTheyIntersect(const LineSegment &m) const {   Point p[4];   p[0] = begin();   p[1] = m.begin();   p[2] = end();   p[3] = m.end();   unsigned char flag = 0;   {          float z  = (p[1].x - p[0].x) * (p[2].y - p[1].y) -                 (p[1].y - p[0].y) * (p[2].x - p[1].x);          if (z > 0)         flag = 2;          else if (z < 0)         flag = 1;   }   {          float z  = (p[2].x - p[1].x) * (p[3].y - p[2].y) -                 (p[2].y - p[1].y) * (p[3].x - p[2].x);          if (z > 0)         flag |= 2;          else if (z < 0)         flag |= 1;          if (flag == 3)         return false;   }   {          float z  = (p[3].x - p[2].x) * (p[0].y - p[3].y) -                 (p[3].y - p[2].y) * (p[0].x - p[3].x);          if (z > 0)         flag |= 2;          else if (z < 0)         flag |= 1;          if (flag == 3)         return false;   }   {          float z  = (p[0].x - p[3].x) * (p[1].y - p[0].y) -                 (p[0].y - p[3].y) * (p[1].x - p[0].x);          if (z > 0)         flag |= 2;          else if (z < 0)         flag |= 1;   }   return flag != 3;    } inline void set(const Point& begin, const Point& end) {        begin_ = begin;        vector_ = end - begin; } inline void setOriginAndVector(const Point& begin, const Point& v) {        begin_ = begin;        vector_ = v; } /* Given this Line, starting from begin_ and moving towards end, then turning towards P2, is the turn clockwise, counterclockwise, or straight? Note: for a counterclockwise polygon of which this segment is a side, Clockwise means P2 would "light the outer side" and Counterclockwise means P2 would "light the inner side". Straight means colinear. */    inline Clockness clockness(const Point& P2) const    {        const float a = vector_.x * (P2.y - begin_.y) - (P2.x - begin_.x) * vector_.y;        if (a > 0)       return Counterclockwise; // aka left        if (a < 0)       return Clockwise;        return Straight;    }    inline bool clocknessIsClockwise(const Point& P2) const    {        return vector_.x * (P2.y - begin_.y) - (P2.x - begin_.x) * vector_.y < 0;    }    //relative to origin    inline Clockness myClockness() const {return clocknessOrigin(begin(), end());}    inline bool clockOK() const {return myClockness() == Counterclockwise;}    //is clockOK(), this is true if p and center are on opposide sides of me //if p is on the line, this returns false    inline bool outside(const Point p) const {   return clockness(p) == Clockwise;    }    inline bool outsideOrColinear(const Point p) const {   return clockness(p) != Counterclockwise;    } void print() const      {        begin().print(); printf(" to "); end().print(); } }; class Wall; /* A wedge is a line segment that denotes a wall, and two rays from the center, that denote the relevant left and right bound that matter when looking at this wall. Initially, the left and right bound are set from the wall's endpoints, as those are the edges of the shadow. But walls in front of (eg, centerward) of this wall might occlude the end points, and we detect that when we add this wedge. If it happens, we use the occluding wall's endpoints to nudge our own shadow rays. The idea is to  minimise the shadow bounds of any given wall by cutting away areas that are already  occluded by closer walls. That way, a given point to test can often avoid being tested against multiple, overlapping areas. More important, if we nudge the effective left and right rays for this wall until they meet or pass each other, that means this wall is completely occluded, and we can discard it entirely, which is the holy grail of this code. Fewer walls means faster code. For any point that's between the effective left and right rays of a given wall, the next question is if it's behind the wall. If it is, it's definitively occluded and we don't need to test it any more. Otherwise, on to the next wall. */ class AggregateWedge; enum VectorComparison {ColinearWithVector, RightOfVector, OppositeVector, LeftOfVector}; static VectorComparison compareVectors(const Point& reference, const Point& point) {    switch (clocknessOrigin(reference, point))    {    case Clockwise:   return RightOfVector;    case Counterclockwise:   return LeftOfVector;    }    if (reference.dot(point) > 0)   return ColinearWithVector;    return OppositeVector; } class LittleTree { public:    enum WhichVec {Left2, Right1, Right2} whichVector; //(sort tie-breaker), right must be greater than left    const Point*   position;  //vector end    LittleTree*        greater; //that is, further around to the right    LittleTree*        lesser; //that is, less far around to the right    LittleTree() {greater = lesser = NULL;}        void readTree(WhichVec* at, int* ip)    {   if (lesser)          lesser->readTree(at, ip);   at[*ip] = whichVector;   ++*ip;   if (greater)          greater->readTree(at, ip);    }    //walk the tree in order, filling an array    void readTree(WhichVec* at)    {   int i = 0;   readTree(at, &i);    } }; class VectorPair { public:    Point    leftVector;    Point    rightVector;    bool acute;    VectorPair() {}        VectorPair(const Point& left, const Point& right)    {   leftVector = left;   rightVector = right;   acute = true;    }    bool isAllEncompassing() const {return leftVector.x == 0 && leftVector.y == 0;}    void set(const Point& left, const Point& right)    {   leftVector = left;   rightVector = right;   acute = clocknessOrigin(leftVector, rightVector) == Clockwise;    }    void setKnownAcute(const Point& left, const Point& right)    {   leftVector = left;   rightVector = right;   acute = true;    }    void setAllEncompassing()    {   acute = false;   leftVector = rightVector = Point(0,0);    }    bool isIn(const Point p) const    {   if (acute)          return clocknessOrigin( leftVector, p) != Counterclockwise &&               clocknessOrigin(rightVector, p) != Clockwise;   //this accepts all points if leftVector == 0,0   return clocknessOrigin( leftVector, p) != Counterclockwise ||              clocknessOrigin(rightVector, p) != Clockwise;    }    //true if we adopted the pair into ourselves. False if disjoint.    bool update(const VectorPair& v)    {   /* I might completely enclose him - that means no change   I might be completely disjoint - that means no change, but work elsewhere   He might enclose all of me - I take on his bounds   We could overlap; I take on some of his bounds   --   We figure this by starting at L1 and moving clockwise, hitting (in some   order) R2, L2 and R1. Those 3 can appear in any order as we move   clockwise, and some can be colinear (in which case, we pretend a convenient order).   Where L1 and R1 are the bounds we want to update, we have 6 cases:   L1 L2 R1 R2 - new bounds are L1 R2 (ie, update our R)   L1 L2 R2 R1 - no change, L1 R1 already encloses L2 R2   L1 R1 L2 R2 - the pairs are disjoint, no change, but a new pair has to be managed   L1 R1 R2 L2 - new bounds are L2 R2; it swallowed us (update both)   L1 R2 L2 R1 - all encompassing; set bounds both to 0,0   L1 R2 R1 L2 - new bounds are L2 R1 (ie, update our L)       If any two rays are colinear, sort them so that left comes first, then right.   If 2 lefts or 2 rights, order doesn't really matter. The left/right case does because   we want L1 R1 L2 R2, where R1=L2, to be processed as L1 L2 R1 R2 (update R, not disjoint)   */   //special cases - if we already have the whole circle, update doesn't do anything   if (isAllEncompassing())          return true; //v is part of this wedge (everything is)   //if we're being updated by a full circle...   if (v.isAllEncompassing())   {          setAllEncompassing(); //we become one          return true;   }   /*   Now we just need to identify which order the 3 other lines are in, relative to L1.        Not so easy since we don't want to resort to arctan or anything else that risks   any roundoff. But clockness from L1 puts them either Clockwise (sooner), or   Straight (use dot product to see if same as L1 or after Clockwise), or   CounterClockwise (later). Within that, we can use clockness between points to sort between them.   */   //get the points R1, L2 and R2 listed so we can sort them by how far around to the right   // they are from L1   LittleTree list[3];   //order we add them in here doesn't matter   list[0].whichVector = LittleTree::Right1;   list[0].position = &this->rightVector;   list[1].whichVector = LittleTree::Left2;   list[1].position = &v.leftVector;             list[2].whichVector = LittleTree::Right2;   list[2].position = &v.rightVector;   //[0] will be top of tree; add in 1 & 2 under it somewhere   for (int i = 1; i < 3; ++i)   {          LittleTree* at = &list[0];          do {         bool IisGreater = list[i].whichVector > at->whichVector; //default if nothing else works         VectorComparison L1ToAt = compareVectors(leftVector, *at->position);         VectorComparison L1ToI = compareVectors(leftVector, *list[i].position);         if (L1ToI < L1ToAt)                IisGreater = false;         else if (L1ToI > L1ToAt)                IisGreater = true;         else         {                if (L1ToI != OppositeVector && L1ToI != ColinearWithVector)                {               //they are in the same general half circle, so this works               switch (clocknessOrigin(*at->position, *list[i].position))               {               case Clockwise: IisGreater = true; break;               case Counterclockwise: IisGreater = false; break;               }                }         }         //now we know where [i] goes (unless something else is there)         if (IisGreater)         {                if (at->greater == NULL)                {               at->greater = &list[i];               break; //done searching for [I]'s place                }                at = at->greater;                continue;         }         if (at->lesser == NULL)         {                at->lesser = &list[i];                break; //done searching for [I]'s place         }         at = at->lesser;         continue;          } while (true);   }   //we have a tree with proper order. Read out the vector ids   LittleTree::WhichVec sortedList[3];   list[0].readTree(sortedList);   unsigned int caseId = (sortedList[0] << 2) | sortedList[1]; //form ids into a key. Two is enough to be unique   switch (caseId)   {   case (LittleTree::Left2 << 2) | LittleTree::Right2: //L1 L2 R2 R1          return true; //no change, we just adopt it   case (LittleTree::Right1 << 2) | LittleTree::Left2: //L1 R1 L2 R2          return false; //disjoint!   case (LittleTree::Right1 << 2) | LittleTree::Right2: //L1 R1 R2 L2          *this = v;          return true; //we take on his bounds   case (LittleTree::Right2 << 2) | LittleTree::Left2: //L1 R2 L2 R1          setAllEncompassing();          return true; //now we have everything   case (LittleTree::Left2 << 2) | LittleTree::Right1: //L1 L2 R1 R2          rightVector = v.rightVector;          break;   default: //(LittleTree::Right2 << 2) | LittleTree::Right1: //L1 R2 R1 L2          leftVector = v.leftVector;          break;   }   //we need to fix acute   acute = clocknessOrigin(leftVector, rightVector) == Clockwise;   return true;    } }; class Wedge { public:    //all points relative to center    LineSegment wall; //begin is the clockwise, right hand direction    Point    leftSideVector;  //ray from center to this defines left or "end" side    Wedge*   leftSidePoker; //if I'm updated, who did it    Point    rightSideVector;  //ray from center to this defines left or "end" side    Wedge*   rightSidePoker;  //if I'm updated, who did it    Wall*    source;        //original Wall of .wall    VectorPair  outVectors;    float   nearestDistance;  //how close wall gets to origin (squared)    AggregateWedge* myAggregate;  //what am I part of?    inline Wedge(): source(NULL), leftSidePoker(NULL), rightSidePoker(NULL), myAggregate(NULL) {}    void setInitialVectors()    {   leftSidePoker = rightSidePoker = NULL;   rightSideVector = wall.begin();   leftSideVector = wall.end();   outVectors.setKnownAcute(wall.end(), wall.begin());    }    inline bool testOccluded(const Point p, const float distSq) const //relative to center    {   if (distSq < nearestDistance)          return false; //it cannot   if (clocknessOriginIsCounterClockwise(leftSideVector, p))          return false; //not mine   if (clocknessOriginIsClockwise(rightSideVector, p))          return false; //not mine   return wall.outside(p); //on the outside }    inline bool testOccludedOuter(const Point p, const float distSq) const //relative to center    {   if (distSq < nearestDistance) //this helps a surprising amount in at least Enya          return false; //it cannot   return wall.outside(p) && outVectors.isIn(p); //on the outside }    inline bool nudgeLeftVector(Wedge* wedge)    {   /*   So. wedge occludes at least part of my wall, on the left side.   It might actually be the case of an adjacent wall to my left. If so,   my end() is his begin(). And if so, I can change HIS rightSideVectorOut   to my right (begin) point, assuming my begin point is forward of (or on)   his wall. That means he can help kill other walls better.   */   if (wedge->wall.begin() == wall.end() && !wedge->wall.outside(wall.begin())) //is it legal?   {          outVectors.update(VectorPair(wedge->wall.end(), wedge->wall.begin()));          wedge->outVectors.update(VectorPair(wall.end(), wall.begin()));   }   //turning this on drives the final wedge down, but not very much   bool okToDoOut = true;   bool improved = false;   do {          if (wall.outside(wedge->wall.begin()))         break; //illegal move, stop here          if (clocknessOrigin(leftSideVector, wedge->wall.begin()) == Clockwise)          {         leftSideVector = wedge->wall.begin();         leftSidePoker = wedge;         improved = true;          }          if (okToDoOut)          {         okToDoOut = !wall.outside(wedge->wall.end());         if (okToDoOut)                outVectors.update(VectorPair(wedge->wall.end(), wall.begin()));          }          wedge = wedge->rightSidePoker;   } while (wedge);   return improved;    }        inline bool nudgeRightVector(Wedge* wedge)    {   /*   So. wedge occludes at least part of my wall, on the right side.   It might actually be the case of an adjacent wall to my right. If so,   my begin() is his end(). And if so, I can change HIS leftSideVectorOut   to my left (end() point, assuming my begin point is forward of (or on)   his wall. That means he can help kill other walls better.   */   if (wedge->wall.end() == wall.begin() && !wedge->wall.outside(wall.end())) //is it legal?   {          outVectors.update(VectorPair(wedge->wall.end(), wedge->wall.begin()));          wedge->outVectors.update(VectorPair(wall.end(), wall.begin()));   }   //turning this on drives the final wedge count down, but not very much   bool okToDoOut = true;   bool improved = false;   do {          if (wall.outside(wedge->wall.end()))         return improved; //illegal move          if (clocknessOrigin(rightSideVector, wedge->wall.end()) == Counterclockwise)          {         rightSideVector = wedge->wall.end();         rightSidePoker = wedge;         improved = true;          }          if (okToDoOut)          {         okToDoOut = !wall.outside(wedge->wall.begin());         if (okToDoOut)                outVectors.update(VectorPair(wall.end(), wedge->wall.begin()));          }          wedge = wedge->leftSidePoker;   } while (wedge);   return improved;    } }; class AggregateWedge { public:    VectorPair  vectors;    AggregateWedge*   nowOwnedBy;    bool       dead;    AggregateWedge() : nowOwnedBy(NULL), dead(false) {}    bool isIn(const Point& p) const    {   return vectors.isIn(p);    }    bool isAllEncompassing() const {return vectors.leftVector.x == 0 && vectors.leftVector.y == 0;}    void init(Wedge* w)    {   vectors.setKnownAcute(w->leftSideVector, w->rightSideVector);   w->myAggregate = this;   nowOwnedBy = NULL;   dead = false;    }    //true if it caused a merge    bool testAndAdd(Wedge* w)    {   if (dead) //was I redirected?          return false; //then I don't do anything   if (!vectors.update(VectorPair(w->wall.end(), w->wall.begin())))          return false; //disjoint   AggregateWedge* previousAggregate = w->myAggregate;   w->myAggregate = this; //now I belong to this   if (previousAggregate != NULL) //then it's a merge   {          vectors.update(previousAggregate->vectors);          //That means we have to redirect that to this          assert(previousAggregate->nowOwnedBy == NULL);          previousAggregate->nowOwnedBy = this;          previousAggregate->dead = true;          return true;   }   return false;    } }; class AggregateWedgeSet { public:    int           at;    int           firstValid;    AggregateWedge    agList[8192];    float         minDistanceSq;    float         maxDistanceSq;    AggregateWedgeSet() : minDistanceSq(0), maxDistanceSq(FLT_MAX)  {}    void add(int numberWedges, Wedge* wedgeList)    {   at = 0;   for (int j = 0; j < numberWedges; ++j)   {          Wedge* w = wedgeList + j;          w->myAggregate = NULL; //none yet          bool mergesHappened = false;          for (int i = 0; i < at; ++i)         mergesHappened |= agList[i].testAndAdd(w);          if (mergesHappened)          {         //some number of aggregates got merged into w->myAggregate         //We need to do fixups on the wedges' pointers         for (int k = 0; k < j; ++k)         {                AggregateWedge* in = wedgeList[k].myAggregate;                if (in->nowOwnedBy) //do you need an update?                {               in = in->nowOwnedBy;               while (in->nowOwnedBy) //any more?                      in = in->nowOwnedBy;               wedgeList[k].myAggregate = in;                }         }         for (int k = 0; k < at; ++k)                agList[k].nowOwnedBy = NULL;          }          if (w->myAggregate == NULL) //time to start a new one          {         agList[at++].init(w);          }   } // all wedges in   minDistanceSq = FLT_MAX;   for (int j = 0; j < numberWedges; ++j)   {          //get nearest approach          float ds = wedgeList[j].nearestDistance;          if (ds < minDistanceSq)         minDistanceSq = ds;   }   minDistanceSq -= 0.25f; //fear roundoff - pull this is a little   firstValid = 0;   for (int i = 0; i < at; ++i)          if (!agList[i].dead)          {         firstValid = i; #if 0 // Not sure this is working? Maybe relates to using L to change bounds?         //if this is the only valid wedge and it is all-encompassing, then we can         //walk all the wedges and find the furthest away point (which will be some         //wall endpoint). Anything beyond that cannot be in bounds.         if (agList[i].isAllEncompassing())         {                maxDistanceSq = 0;                for (int j = 0; j < numberWedges; ++j)                {               float ds = wedgeList[j].wall.begin().dotSelf();               if (ds > maxDistanceSq)                      maxDistanceSq = ds;               ds = wedgeList[j].wall.end().dotSelf();               if (ds > maxDistanceSq)                      maxDistanceSq = ds;                }         } #endif         break;          }    }    const AggregateWedge* whichAggregateWedge(const Point p) const    {   for (int i = firstValid; i < at; ++i)   {          if (agList[i].dead)         continue;          if (agList[i].isIn(p))          {         return agList + i;          }   }   return NULL;    } }; //#define UsingOuter //this slows us down. Do not use. #ifdef UsingOuter #define TheTest testOccludedOuter #else #define TheTest testOccluded #endif class AreaOfView { public:    Point    center;    float    radiusSquared;    int      numberWedges;    BoundingRect  bounds;    Wedge    wedges[8192];    //VERY experimental    AggregateWedgeSet ags;    inline AreaOfView(const Point& center_, const float radius) :        center(center_),        radiusSquared(radius * radius),        numberWedges(0)      {        bounds.set(center, radius);        addWalls(); }     void changeTo(const Point& center_, const float radius)    {   center = center_;        radiusSquared = radius * radius;   bounds.set(center, radius);   numberWedges = 0;   addWalls();    }    void recompute() //rebuild the wedges, with existing center and radius    {   bounds.set(center, sqrtf(radiusSquared));   numberWedges = 0;   addWalls();    }    inline bool isIn(Point p) const    {   p -= center;   const float distSq = p.dotSelf();   if (distSq >= radiusSquared)          return false;   for (int i = 0; i < numberWedges; ++i)   {          if (wedges[i].TheTest(p, distSq))         return false;   }   return true;    }    /* On the theory that the wedge that rejected your last point has a higher than    average chance of rejecting your next one, let the calling thread provide    space to maintain the index of the last hit    */    inline bool isInWithCheat(Point p, int* hack) const    {   p -= center;   const float distSq = p.dotSelf();   if (distSq >= radiusSquared)          return false;   if (distSq < ags.minDistanceSq)          return true; //this range is always unencumbered by walls   if (distSq > ags.maxDistanceSq) //not working. Why?          return false;   if (numberWedges == 0)          return true; //no boundaries   //try whatever worked last time, first. It will tend to win again   if (wedges[*hack].TheTest(p, distSq))   {          return false;   } #define UseAgg #define UseAggP #ifdef UseAgg   const AggregateWedge* whichHasMe = ags.whichAggregateWedge(p);   if (whichHasMe == NULL)          return true; //can't be occluded! #endif   //try everything else   for (int i = 0; i < *hack; ++i)   { #ifdef UseAggP #ifdef UseAgg          if (wedges[i].myAggregate != whichHasMe)         continue; #endif #endif          if (wedges[i].TheTest(p, distSq))          {         *hack = i; //remember what worked for next time         return false;          }   }   for (int i = *hack + 1; i < numberWedges ; ++i)   { #ifdef UseAggP #ifdef UseAgg //does seem to help speed, but don't work yet          if (wedges[i].myAggregate != whichHasMe)         continue; #endif #endif          if (wedges[i].TheTest(p, distSq))          {         *hack = i; //remember what worked for next time         return false;          }   }   return true;    }    inline bool isInWithWallExclusion(Point p, const Wall* excludeWall) const    {   p -= center;   const float distSq = p.dotSelf();   if (distSq >= radiusSquared)          return false;   for (int i = 0; i < numberWedges; ++i)   {          if (wedges[i].source == excludeWall)//this one doesn't count         continue;          if (wedges[i].TheTest(p, distSq ))         return false;   }   return true;    }    void addWall(Wall* w, const float nearestDistance);    void addWalls(); }; class AreaRef { public:    AreaOfView*    a;    AreaRef() {a = NULL;}    void set(const Point& p, float radius)    {   if (a == NULL)          a = new AreaOfView(p, radius);   else          a->changeTo(p, radius);    }    ~AreaRef() {delete a;}    void empty() {delete a; a = NULL;}    AreaOfView* operator->() const {return a;} }; class WallSet { public:    int   length;    int   at;    WallAndDist*   list;    WallSet() {   at = 0;   length = 2038;   list = (WallAndDist*)malloc(length * sizeof(*list));    }    ~WallSet() {free(list);}    void add(Wall* w, const float distSq)    {   if (at >= length)   {          length *= 2;          list = (WallAndDist*)realloc(list, length * sizeof(*list));   }   list[at].wall = w;   const LineSeg* s = w->getSeg();   list[at].lenSq = s->p[0].distanceSq(s->p[1]);   list[at++].distSq = distSq;    }    inline void sortByCloseness()    {   qsort(list, at, sizeof *list, cmpWallDist);    } }; void AreaOfView::addWall(Wall* w, const float nearestDistance) {    if (num

Cellular Textures, the light speed approach

Mon, 10 Aug 2009 16:34:36 +0000

Introduction
Hi, my name is Carsten Przyluczky (aka Kastor formerly known as Jazzoid), I am a student of computer science. I am interested in game development and a huge fan of the demo scene, every time I saw a good Demo in 64 K or even 4 K I asked myself: “How the heck do they do that?! How can they put that much content into such a small space?” So me and a friend came up with the idea of writing a complete game in 96 K. During my research, I found out that most of the content is “procedural content”. That means it is calculated at runtime. Yes, they do some Math-magic and the results are these good-looking textures, music or even 3D models.

Creating procedural textures is mostly done by taking some basic shapes or patterns, and manipulating them. Manipulations can be colour inversion, blending or adding to textures, and so on. Basic shapes are things like circles, rectangles, plasma noise, or cellular patterns / textures.

In this article I am going to explain how to generate these cellular textures. For those who have never seen a cellular texture they look like this:

They are very versatile, you might use them to create spider webs or cracked dirt textures. If you search the internet you will find some great tutorials (see link section). They mostly involve an easy but slow approach. I will explain this basic approach too, and improve it step by step. In the end we will have an algorithm that is much faster and has some nice extra features.

The basic approach (turtle speed)
Well, the basic idea is easy. Let's assume you want to create a texture that is 256 x 256 pixels large. The first step of the algorithm is to generate a set of random points let's call this ^P. As you can see, the x and y coordinates must lie between 0 and texture size – 1. So, for each point ^q of the texture you have to find the point p∈^P that has the smallest distance to ^q, and store it into an array. In the case of q∈^P the distance is set to zero.

The distance we just found will determine the color of the pixel in our texture. Usually, the distance is calculated by using the Euclidean formula. For our example it would look like this:

There are some more possibilities, but let's skip these for now. Once you have calculated all distances, you have to normalise the values in the array to the range of 0 to 255.

These values can be easily transformed into a grey scale value by setting red green and blue to the value of the normalized array.

That's pretty simple, right? Just calculate the distances and use them to colourize your texture.

As far as I know this is the most common way to create the nice-looking cellular textures. Well, the main part of our algorithm is to find the closest point of our random set. If we want to speed up the algorithm we need to speed up that part.

But first, it's time to have a look at the algorithm's running time.

Let ^m be the cardinality of ^P, our random point set, and let be ⁿ the count of points in our texture. It is clear that the common algorithm needs to calculate distances for each of our points. This leads to a running time of,

O(n * m)

This means, if you want your texture to look finer and you use more random points, the time it takes to calculate rises drastically.

The BSP approach (cheetah speed)
The first optimization that came into my mind was to use a binary space tree, to store the random points. Once we've created our tree we can use it to speed up “the quest for the closest point” In a BSP, each inner node holds two points and each leave one. In addition, each inner node has two branches, so we need at least two points to build a tree. The left branch of an inner node holds points that are closer to the first Node. The right branch works the other way around. So, if you search the closest point of our random set, you start at the root, calculate which point is closer and follow the corresponding branch.

Ok, lots of information here, let's take an example:

The red points are our random points. The first point added is A then B , C and so on. If we add the points in another order the tree will look different. You need to check if the tree is balanced, in the worst case you will deal with a regular list, which won't be faster.

This approach increases the speed but needs some preprocessing like rebalancing and so on. For me this is no option because it needs too much extra code. Nevertheless, the running time of this approach decreases to something like

O(n * log(m))

(average case).

The grid approach (light speed)
OK, fasten your seatbelts - we are reaching light speed now!

First of all remember that the color of a texture point flows from the distance to the closest point of our random set ^P. So, that's our main task. If we want to speed up our algorithm, optimizing this task is a good start.

We can archive this using a grid. We put the grid over our texture. Let's start with a 8 x 8 grid:

As you can see our 8x8 grid has 64 cells. Actually the grid is just a virtual construct to relax the problem. The main trick now is to pick a random point for each cell. So in our case we need 64 random points. The x and y values of our points must lie in the range of the corresponding cell boundaries, like this

The blue points represent our random set ^P. Okay now look at this:

Let the red point be an arbitrary point q, and we want to know the closest point p∈^P. All we need to do is calculate the distances to the point that shares the same cell as q, and the eight surrounding ones. In the image the corresponding cells are grey. That's it ! One of these eight points must be the closest. Any other point must have a greater distance.

After only nine distance calculations, we got our closest point p∈^P! Now that the cardinality of ^P doesn't matter any more, the running time of this approach is,

O(n * 9)

no matter whether you take 16 or 1024 random points!

Extra features

Making it tillable
Often you want your textures to be tillable. To achieve this you need to consider the “wrapped” points too. Let's call the cell that contains our arbitrary point ^q Q-Cell. And now let's assume Q-Cell is the upper left cell. The Q-Cell has no left neighbor, so we need to wrap, and take the upper right Cell as the Q-Cell's left neighbor. Look at the following image, it should make things clear: In my code I do something like this to perform the wrap (Pseudo code):

LeftNeighborIndex = (QCell.X+NumOfCellsPerRow-1) mod NumOfCellsPerRow I perform a kind of “shift” here, and the -1 should “give” me the left-hand neighbor. In our case the Q-Cell has the (horizontal) index 0 and its right-hand neighbor (far-right) the index 7. If we put this information into our formula the result is

(0+8-1) mod 8 = 7 → correct! And as another example we assume our Q-Cell is the upper right-hand, and we want to know what its right-hand neighbor is (add 1) our calculation looks like this

(7+8+1) mod 8 = 16 mo 8 = 0 → correct!
Create different shapes
If you look at the three images in the introduction, you will notice that the left-hand one looks different. I used another distance function here, namely the distance of the two closest points. As I mentioned before, you can use several distance functions. Play around with it, invert it, square it, multiply the two closest, there are no limits! The following two images are created with Manhattan-noise. It is realized by taking the two closest points and subtract their x and y values

Flaws
Well like all good things, this approach has some flaws. There are some extremly rare cases when this algorithm produces not the same results as the naive one. Take a look at this picture:

  The green line indicates the smallest distance. As you can see the grid algorithm will select another point. This can be fixed by widen the grid and check 5x5 instead of 3x3. I can tell you that I did a lot of tests with 3x3 grid and the results were equal, but I still want to point that weakness out.

Closing Words
So, this is the end my friend. I hope you had fun reading this article, because it was a lot of fun writing it. Go to my blog to find out more about my work and my procedural content editor (see link section). If you have any questions feel free to write an email to jazzoid@gmx.de. Cheers,
Carsten

Links
http://kastor.wordpress.com/
http://graphicdesignertoolbox.com/screenshots (look at worley noise)
http://www.blackpawn.com/texts/cellular/default.html
http://petewarden.com/notes/archives/2005/05/testing.html
http://www.worley.com/

Image Space Lighting

Thu, 11 Jun 2009 14:23:51 +0000

Introduction
Current graphics hardware has built in functionality for rendering a constant number of light sources at a time. While for some scenes this constant is large enough, for many scenes it is not, and for those scenes this hardware-determined limit can be frustrating to work around. To increase the number of light sources, most developers end up having to use techniques that only approximate the lighting (like glow effects) or are slow (like looping through all the lights manually per pixel or per vertex). This article focuses on an efficient way to render hundreds of dim light sources accurately using an image space technique that takes advantage of the graphics processor’s fast rasterization. This technique is very useful for scenes that have a lot of lampposts, lanterns, towns at night, lights on models, light-emitting particle effects and other scenes that have a large number of ‘small’ light sources. Additionally, since this technique lies in image space, the cost of the algorithm does not depend on the complexity of the geometry in the scene, only the light sources themselves. The same light sources rendered with this technique have the same additional cost whether the scene has 10K or 100K polygons.

In order to run this algorithm efficiently, the graphics card must support multiple render targets. This means that this algorithm only works for graphics cards that support DirectX 9.0 or OpenGL 2.0 (you can use extensions for MRT in earlier versions). Before reading this article, you should understand how to render to multiple targets and how to program vertex and pixel shaders.

Lighting
Before we jump into the algorithm itself, lets first briefly explain lighting, and in doing so develop intuition on how the algorithm can work efficiently. For this algorithm, we only consider lights that actually have a position and exist in world space (unlike directional lights). To render directional lights, the built-in light sources on the graphics card should suffice (as there usually are not too many directional lights). For simplicity, we will assume all light sources are point light sources, although there are simple extensions to the algorithm to render non-point light sources which I will briefly discuss later.

Let’s say there are three point lights in a scene.


Figure 1: Three light sources in a simple scene. The black suns represent the actual light sources, and the white circles around them represent their area of effect.

  These light sources cast light in all directions around them, however the power of the light source attenuates quadratically with distance. This means that the farther away from the light source a point is, the less light it receives. At some distance, the amount of light received from the light source is small enough (less than some delta) that its contribution can be considered zero. This distance defines a sphere around the light source that represents the light source’s area of effect. Everything within this sphere receives light from the light source, and everything outside this sphere receives almost no light from the light source.

When rendering light sources using the hardware’s built in light sources or manually looping through a list of light sources, this area of effect is ignored, and all pixels will calculate the contribution of every light source, no matter how far away it is. If the light source is large and its area of effect covers the entire screen, then no computation is wasted this way. But when rendering thousands of dim light sources (which could have an area of effect of only a few pixels), this approach ends up wasting computational power for all of the pixels where the light’s contribution is negligible (outside of the light’s area of effect).

If you have 1,000 point lights, each of which affects only 100 pixels on the screen (where a screen might have 1,000,000 pixels), then without considering the area of effect, you are computing lighting for 1,000 lights over 1,000,000 pixels, a total of 1,000,000,000 lighting calculations. However, if for every light source, you only compute the lighting for the pixels within that light’s area of effect, then you only compute lighting for 100 pixels for each light source, for a total of 100,000 lighting calculations. This is a 10,000x speedup. Obviously, taking the light source’s area of effect into account is a big win when rendering lots of small light sources – and that is the essential idea behind image space lighting, and how it can operate so efficiently.

The Area of Effect
Calculating the area of effect of a point light is a simple computation. The contribution of a light source on graphics hardware can be defined by the following:  Contribution = Light Power * Cosine Term / (d ^2 * QUADRATIC_ATTENUATION + d * LINEAR_ATTENUATION + CONSTANT_ATTENUATION)

  Where the light power is the power (or brightness) of the light source (a constant for each light source), The cosine term is unknown, as it depends on the normal of the point being lit, but it is strictly between zero and one, so we conservatively set it to one. The three attenuation parameters define how the light’s power attenuates over distance (where d is the distance between the light source and the point). These are constants for each light source. We can then rearrange this equation to get the following:

  d ^2 * QUADRATIC_ATTENUATION + d * LINEAR_ATTENUATION + CONSTANT_ATTENUATION = Light Power / Contribution

  The equation is now quadratic with respect to d. To solve this equation, we set contribution variable to the smallest value we want to include in our area of effect, and solve for d. The positive solution of this equation is now the distance at which the contribution is equal to our small delta value, and beyond that distance, the contribution decreases (attenuates), so this distance defines the area of effect of the light source. We can compute this distance once for each light source (and update it if the light source changes). Using this distance, we can quickly determine if any point is within the area of effect of the light source or not (by comparing this distance with the distance from the point to the light source).

Proxy Shapes
We can now loop through all light sources per pixel, and test if that pixel is within the light’s area of effect before computing the lighting. However, this will not greatly increase the speed – as branching on the GPU is slow, and a loop with an if statement for 1000 light sources for every pixel is still very expensive. Intuitively, we want to quickly figure out which pixels in image space are within our light source’s area of effect, and then add the contribution of the light source only for those pixels. Luckily – there is a very fast and optimized algorithm for projecting shapes into screen space on the GPU: Rasterization. GPU rasterization can rasterize over a million triangles in real time. If we bound the light’s area of effect within a low poly object – say a cube (12 triangles), we can then rasterize a cube for each light source to the screen, and the pixels in that cube represent approximately the pixels within the light source’s area of effect (if it is a bounding shape, it will never under-approximate).


Figure 2: Bounding cubes rasterized for each point light, quickly generating pixels within the light’s area of effect. Proxy cubes shown as transparent boxes, transparent circles depict the actual area of effect projected on the ground.

  If the GPU can render a million triangles and each light source’s proxy shape is a cube (12 triangles), then we can potentially render over 80,000 light sources! Practically, we cannot achieve anywhere near 80,000 light sources, but this shows how using rasterization can easily and efficiently generate the pixels within each light source’s area of effect.

Rendering and Implementation
Now that we know the theory behind image space lighting, it is time to work out the actual rendering and implementation details. The steps to rendering image space lighting are as follows:
Extract the position, normal and color information of the scene to calculate the lighting at any pixel. This can be done with deferred shading, where the scene is rendered once with a pixel shader that stores the position, normal and color information into the proper render targets. Any effects like normal mapping can be done here – just outputting the final normal to the normal texture.
Render the scene to the screen normally (as it was rendered off screen previously). However, this does not require rendering the scene geometry again, as the per pixel attribute results were already generated in the deferred shading pass. Simply render a quad that fills the entire screen (this is easy to do with a vertex shader), and then retrieve the per pixel position, normal and color information from the deferred shading textures. During this pass, any primary light sources (like directional lights) and effects can be rendered (the contribution to the scene from the point lights will be added later). It is worth noting that the pixel shader only executes once per pixel in this pass, which can improve the speed of the rendering for complex pixel shaders.
Render the point lights on top of the previously rendered scene as follows:

Bind a shader that computes the lighting contribution from the light source whose proxy shape is being rendered on the point at the rasterized pixel position in the deferred shading textures. This shader can also discard the fragment if the point is not within the light’s area of effect to save computation costs (which happens because the proxy shape over-bounds the area of effect and because the depth test is disabled when rendering the proxy shapes (see below), so the point could be far away in the depth coordinate).
Set the blending function to 1 * src + 1 * dst. This indicates that the computed lighting at each pixel is added to the previous value. As the lighting equation is additive (summation of the contribution of all light sources), this works as expected.
Enable face culling, and cull front faces. If the back faces are culled, nothing will be rendered when the camera moves inside the proxy shape (as the front faces will be behind the camera). As long as the proxy shape is convex, there are only two rasterized fragments for each pixel on the shape (one on the front face, and one on the back face). If the front faces are culled, only one fragment is rasterized per pixel, which means that no light will contribute lighting twice for any pixel (which would be incorrect).
Disable the depth test as pixels on the back face of the proxy shape could fail the depth test (thus contributing no light) even when the light source should contribute light to that pixel. However, since the pixel shader can quickly discard pixels outside of the light source’s area of effect, this does not affect performance too much.
Render all of the proxy shapes. The proxy shapes can be stored in a tree structure to quickly render only the proxy shapes within the view frustum to improve performance.
Finally, unbind the shader and textures, and then set blending, face culling and depth test back to their previous values.
Code Snippets
Below is an OpenGL code snippet from the game project Aero Empire [2], which renders the point lights following the above implementation.  Shader* cur_shader; void renderLight(Light* light){   //pass light world position to pixel shader   cur_shader->setUniformVec3("lightWorldPos", light->getWorldPosition());   //pass light color (magnitude of vector is power) to pixel shader   cur_shader->setUniformVec3("lightColor", light->getColor());   //pass the light's area of effect radius to pixel shader   cur_shader->setUniformFloat("lightAoE", light->getAoERadius());    light->render(); //render the proxy shape } void renderLights(const Shape* scene){   //set blend function   glBlendFunc(GL_ONE, GL_ONE);   //cull front faces   glEnable(GL_CULL_FACE);    glCullFace(GL_FRONT);      //disable depth testing   glDisable(GL_DEPTH_TEST);    glDepthMask(false);      //bind point light pixel shader   Shader* shader = shaders[SHADER_LIGHT];    shader->bind();   //bind position, normal and color textures from deferred shading pass   bindTexture(positionMap, POSITION_UNIT);    shader->setUniformInt("positionMap", POSITION_UNIT);   bindTexture(normalMap, NORMAL_UNIT);   shader->setUniformInt("normalMap", NORMAL_UNIT);   bindTexture(colorMap, COLOR_UNIT);   shader->setUniformInt("colorMap", COLOR_UNIT);   bindTexture(attrMap, ATTR_UNIT);   shader->setUniformInt("attrMap", ATTR_UNIT);   //set pixel shader attribute "camPosition" to the camera pos    //(as it is needed for Phong shading)   camera.loadCameraPosition(shader, "camPosition");    shader->setUniformFloat("d_sx", 1.0/width);   shader->setUniformFloat("d_sy", 1.0/height);   cur_shader = shader;      //run renderLight function on all light proxy shapes in scene   getLights(scene, &renderLight);   //unbind and reset everything to desired values   shader->unbind();    unbindTexture(POSITION_UNIT);   unbindTexture(NORMAL_UNIT);   unbindTexture(COLOR_UNIT);   unbindTexture(ATTR_UNIT);   glDisable(GL_CULL_FACE);   glEnable(GL_DEPTH_TEST);   glDepthMask(true); }  Below is the GLSL pixel shader bound above which computes the lighting (the vertex shader only sets gl_Position = ftransform() ).  //input parameters uniform sampler2D positionMap, normalMap, colorMap, attrMap; uniform vec3 camPosition, lightWorldPos, lightColor; uniform float lightAoE, d_sx, d_sy; void main(void) {   //calculate screen coord     vec2 coord = vec2(gl_FragCoord.x*d_sx, gl_FragCoord.y*d_sy);               //get the position from deferred shading     vec4 position = texture2D(positionMap, coord);               //vector between light and point     vec3 VP = lightWorldPos-position.xyz;               //get the distance between the light and point     float distance = length(VP);               //if outside of area of effect, discard pixel     if(distance > lightAoE) discard;               //normalize vector between light and point (divide by distance)     VP /= distance;               //get the normal from deferred shading     vec4 normal = texture2D(normalMap, coord);               //get the color from deferred shading     vec4 color = texture2D(colorMap, coord);               //get lighting attributes from deferred shading     vec4 attributes = texture2D(attrMap, coord);          float diff_coefficient = attributes.r;     float phong_coefficient = attributes.g;     float two_sided = attributes.b;     float cos_theda = dot(normal.xyz, VP);                  //calculate two sided lighting.     cos_theda = (cos_theda < 0.0)?-two_sided*cos_theda:cos_theda;               //calculate diffuse shading     float diffuse = diff_coefficient*cos_theda;               //calculate half vector     vec3 H = normalize(VP+normalize(camPosition - position.xyz));               //calculate Phong shading     float phong = phong_coefficient*pow(max(dot(H, normal.xyz), 0.0), 100.0);               //calculate light contribution with attenuation     vec3 C = lightColor*(color.rgb*diffuse+phong)/(distance*distance+0.8);               //all lights have constant quadratic attenuation of 1.0, with a constant attenuation of 0.8 to avoid dividing by small numbers     gl_FragColor = vec4(C, 1.0); //output color }  If you are planning to add lanterns or lamps, adding two-sided lighting allows for point lights inside a lampshade or lantern to contribute light to it. The above shader includes a very simple implementation of two-sided lighting.
Results

Figure 3: 2175 bright cube lanterns rendered on a flat plane in a night scene.
Figure 4: A daytime rendering of a blimp from Aero Empire [2] with lanterns during the day. Notice the Phong lighting on the envelope.

Figure 5: A nighttime rendering of the same blimp as above.

Figure 6: Ship hanger interior from Infinity [3]. This rendering includes both point lights and ambient lights (adding point lights that approximate indirect lighting).

  As you can see, image space lighting creates accurate lighting effects that add to the scene whether daytime or nighttime, interior or exterior. The ability to efficiently render a large number of point light sources allows for much more diverse lighting environments which greatly improves the rendering of the scene, yet still runs in real time.

Performance
The cost of this algorithm depends solely on the number of point lights rendered and the number of pixels within each light’s area of effect. The number of lights determines the overhead of the rasterization step. Without this overhead, in the worst case where the area of effect of all lights is the entire screen, the algorithm would simply add the contribution of every light for every pixel (which would have the same cost as rendering all of those lights by looping over all lights per pixel). However, there is also the overhead of the rasterization step, which is small (as modern rasterizers can handle a million triangles), but is worth noting. The number of pixels rasterized from the proxy shapes determines the main performance hit for this algorithm. Point lights that are very bright, or are very close to the camera end up generating a large number of pixels. While one or two point lights filling the entire screen does not impact performance too much, having many of such point lights is the main bottleneck of image space lighting. This algorithm performs better when the point lights are not clustered, for if the camera is close enough to one point light for it to fill the screen, then the other point lights are farther away and so have a smaller area of effect in screen space (as long as the points are not clustered).

To increase performance out of this algorithm – reduce point light usage, decrease point light brightness, and spread point lights out more (you can always approximate a cluster of point lights with a single brighter point light).

Below are some performance results of this algorithm. All results are from rendering a lantern scene (like Figure 3) at a resolution of 600x400 on a GeForce 9200M GS. The ground is a 2000x2000 unit plane, and lanterns (12 triangles with a point light inside) are scattered in 50x50 unit clusters that are uniformly distributed on the ground plane. The camera is guaranteed to be close to one of those clusters (worst case). The light shader computes two-sided lighting for diffuse surfaces, and the deferred shading pass renders to four 32-bit floating-point textures.

Clusters is the number of clusters generated, lights is the number of lights per cluster generated, ISL # is the fps of the image space lighting algorithm where all lights have an area of effect radius # units (where the area of effect is determined by the brightness of the light and the cutoff threshold), and PPS is the fps of a naive direct loop of all lights per pixel (per pixel sum).



  This shows that if the lights are not very clustered or bright, the image space lighting technique can render hundreds of point lights without too much of a performance hit. In the worst case, where the camera is near a single cluster with many point lights, the image space lighting technique still outperforms the per pixel sum algorithm when the lights are dim, however as the lights get brighter, the image space lighting algorithm slowly converges towards the speed of the per pixel sum algorithm (as it has to compute the lighting for each light source for each pixel).

Downsides and Possible Solutions
There are three major downsides to Deferred Rendering (the first step to Image Space Lighting). All of these downsides come from the fact that you have to store all necessary geometric attributes into textures. The first downside is just the sheer amount of memory storage required (which can be quite limited, depending on the graphics card). In the demo code, I store the position, normal, diffuse color, diffuse component, specular component, and two-sided lighting component. This ends up being 9 32 bit floating point values, 1 32 bit depth value, and 4 8 bit unsigned byte values, for a total of 44 bytes per pixel. For a standard screen size of 1024x768, that is 786,432 pixels and a total of 33 megabytes. When a graphics card only has 128 megabytes dedicated memory, that’s about a fourth of the memory space used, which can be a problem when you also need to store vertex buffer objects (triangle mesh data stored on the GPU) and textures for the various objects in the scene. The storage of geometric data into textures is called a G-Buffer (standing for geometric buffer), and there are several tricks to bring its size down, like using 16 bit floating point values, spherical coordinates, storing the position as just a depth value, etc. In general, you want to compress and pack the data in a way that saves a lot of memory, however, this can also cause banding and other artifacts if you aren’t careful.

The second downside is anti-aliasing. Since you only store the geometric attributes per pixel, you lose the hardware anti-aliasing that allows you to blend edges by computing sub pixels. This causes aliased jagged edges that you have probably seen if you have rendered a scene without anti-aliasing. There are several ways to add anti-aliasing to a deferred rendering approach, like rendering to a larger buffer and then downsampling (however, this requires even more memory and the computational cost of averaging), and doing edge detection and blurring around the edges (which adds computational cost and does not take into account the sub pixels, so it’s not a very effective form of anti-aliasing). The ideal solution would be to generate and store sub pixels only where needed (like at the edges, which wouldn’t increase the memory cost too much), but this would require hardware support and extensions, as graphics card typically only render to a framebuffer (2D texture or buffer).

The final downside is alpha blending. Related to anti-aliasing, sometimes you need to generate multiple samples for certain pixels, but by rendering to a texture, you lose this information (as you can only have one sample per pixel). With alpha blending, you want a semi-transparent object to blend on top of the object behind it, however, this requires the geometric data of the foreground and background object to be stored in the same pixel. Again, with hardware support and extensions for rendering to structures other than framebuffers, you could store all of these samples and easily render scenes with alpha blending. However, without this hardware support, the best way to add alpha blending to deferred rendering scenes is to render the opaque surfaces first, and then render the alpha-blended surfaces on top (however, if they are rendered on top, then they do not get any light contribution from the point lights).

Other Light Types and Effects
Adding spotlights and non-point lights simply requires calculating that light source’s area of effect, and computing a bounding proxy shape for it. For non-point lights, this is a little more complex, as the area of effect may not be spherical. Additionally, the pixel shader would have to handle rendering the different lighting types. Other than that, the algorithm should work the same – for each pixel, add the contribution of all light sources where the point lies within the light source’s area of effect. I have not experimented with non-point lights, but if anyone tries, feel free to share the results. Adding area light sources is also possible by approximating the light with many small point lights (a common approximation for rendering area lights). This allows for light sources that have different shapes (long, square, etc).

Additionally, instead of rendering light sources, you can use this technique to render light patches (for indirect illumination). The C. Dachsbacher and M. Stamminger’s paper, “Splatting Indirect Illumination” [1] uses rasterization to compute surface patches that receive direct lighting, and then computes the lighting received from those patches for all pixels within that patch’s area of effect (as done in image space lighting). This quickly renders an approximation to one bounce indirect illumination.

Currently, there are difficulties in casting shadows from these light sources efficiently. Sampling occlusion from a large number of light sources is a tough and expensive problem – and is usually solved by generating shadow maps for clusters of lights. However, since the lights rendered by this algorithm are dim, the fact that they do not cast shadows is not incredibly noticeable.

Conclusion
Any game that wants to simulate complex lighting effects from a large number of light sources could benefit from this technique. Very few games use more than the hardware’s built-in light sources, and most of the smaller light sources are approximated with glow or light maps. This technique is offered as an alternative to glow and light maps that is still efficient, yet computes accurate lighting effects. This technique would be especially impressive for games with night scenes that are lit up by many small point lights. As a real world example of image space lighting, it is implemented in the game Aero Empire (which I am currently developing). The reason for adding it is mainly for rendering lanterns, which will be common within towns and on the blimps. They will allow for operation of blimps at night, towns lit up by lanterns at night, and lights for the interiors of blimps and buildings (which would be dark otherwise due to shadows). This is one example of how image space lighting can be used in a game. It is my hope that other games will use this technique to add more complex and interesting lighting environments to their graphics programs as well.

Infinity: The Quest for Earth is another real world example of a game that has included image space lighting. F. Brebion, developer of Infinity, implemented a technique similar to image space lighting to improve the interior of hangers and hulls of the ship – which would otherwise be dark when only including light from the sun. The inclusion of indirect lighting also greatly improved the realism of the rendered scenes. See Brebion’s journal on deferred lighting for more information on this [3]. Infinity is another great example of how image space lighting can be used to efficiently add more complex lighting environments.

If you have any questions, improvements or results from using this algorithm that you would like to share, feel free to post in the comments or contact me at terra0nova@hotmail.com.

Further Reading / Sources
I have included the source code for rendering a lantern scene at night with the image space lighting technique. Feel free to fiddle with the code and put together more complex scenes as well. This code should be considered free to use and modify, so long as you give proper credit where due. [1] C. Dachsbacher , M. Stamminger. Splatting Indirect Illumination. 2006.

[2] Collaborative Game Project. Aero Empire. 2009.

[3] F. Brebion. Deferred Lighting and Instant Radiosity. 2009.

Rendering Water as a Post-process Effect

Thu, 04 Jun 2009 19:08:55 +0000
Water plays an important role in all terrain renderers used in modern games and visualizations to present outdoor areas. This is because it improves general image quality in an extent bigger than any other technique and makes it more photorealistic. Alas, even though it tends to look better and better with every game released it is still far from being realistic. Besides, deferred shading gains popularity every month. Many newer games and engines make use of it, eg. Starcraft 2 and Tabula Rasa. Though deferred shading is typically used to limit lighting-related operations from O(objects_number * lights_number) to O(objects_number) I will prove that it can be helpful in many more algorithms.

In this article I will describe a technique of realistic and flexible water rendering by using bump mapping in the post-process stage just after deferred shading. The presented technique fits into the concept of deferred shading pretty nicely, i.e. it helps to avoid additional geometry rendering. A description on deferred shading can be found in [1], [2] and [3]. It is possible to implement this technique using forward rendering but in this case it seems less natural and may require additional work.

The presented algorithm eliminates the majority of the flaws of typical techniques used for water rendering such as hard edges or unrealistic colour extinction. I will talk shortly about common drawbacks in the section “Traditional approaches to water rendering”.

Theory behind water

The theory behind water is very complex and not fully understood. It will suffice to say that there is no adequate formal model describing it so far. Survey conducted by Guillot [4] proved that none of the 46 models he analyzed is valid when compared to reality. Existing models are also a way too complex and computationally expensive to be used for real time applications, especially games. A game should not spend most of its CPU or GPU processing time just to update and render realistic water.

So a completely different solution has to be found.

In my opinion the theory of water for real time applications can be divided into two categories:
Waves - their animation, propagation and interaction with the rest of the world
Optics.
To simulate waves propagation, FFT (Fast Fourier Transform [5]) is frequently used. In this article I do not focus on waves propagation and so I will skip the theory behind it. I strongly believe that what is more important for good looking water is optics and I will describe it in a very detailed way.

Reflection and refraction

A light ray going through the water surface gets reflected and refracted, causing specular, caustics and light shafts to appear. Light shafts and caustics play an important role in underwater scenes, however for an observer standing above the water surface they do not improve the quality of water that much. A key to describing reflection and refraction is the Fresnel term. The Fresnel equation tells us how much light is reflected and refracted when it crosses a border between different media (in this case water and air). We have to make an assumption that light goes through two materials only as the term shown later does not apply in a more general case. Besides, both of them have to be homogeneous and one of them has to have a higher density. "Air – water” pair falls well into this assumption. However "air - window glass” does not as it consists of three materials (air – glass – air). Graphically it is depicted as follows:

Light ray reflection and refraction

For any given angle of light a and refraction coefficients describing both materials A and B, namely n_A and n_B, we use following notation:

The Fresnel term could then be defined as follows:

However this equation is too complex to be computed on GPU per-pixel. Therefore, a frequently used approximation is given below:

This function is very similar to the one described previously so the loss of quality is really insignificant. R(0) is constant and therefore it should be computed only once and then passed to the pixel shader.

Index of refraction for water is equal:

IOR = 1,33333

Colour extinction

One of the most important aspects of proper water rendering are colour extinction as light goes deeper and deeper and also light scattering. However, rarely does anyone pay attention to that, although these are the two phenomena which cause water to have the colour we see. This colour is called the “apparent one” and can be totally different from so-called true water colour. To find this true colour, one has to filter out all the particles and organic material from water. It is done in laboratories and for us is not usable. One thing to note is that true water colour is bluish (not completely transparent as it was believed some time ago). We will use this information later on. For each component of the light spectrum, speed of extinction is different. This is due to different wavelengths of the components. The visible range of light is exceptional in that attenuation fades slowly compared to infra-red or ultraviolet. For them, absorption rate is so high they fade away entirely at a maximum depth of a few centimetres. In the table underneath, average depths at which light components die out are presented. The table is valid only for very clean waters. In the case of muddy ones like lakes or rivers it has to be modified:

We make an assumption that extinction is linear with depth. In general this is not true as water is not homogeneous at its whole depth. Many waters consist of tiers of different temperatures and density which makes them different mediums. However, to simulate water in a convincing way it is not necessary to take this into account.

Besides, extinction speed ratio and water colour are also influenced by chemical composition of the water area bottom, level of siltation, existence of organic materials (like algae or plankton) and even sky colour!

An important conclusion from this deliberation is that the further from the equator the water is the more green its colour becomes.

Traditional approaches to water rendering

The traditional approach to the water rendering problem [6] is based on rendering a plane with reflection and refraction textures applied to it in the level determined by the result of the Fresnel formula. To give water surface the impression of being wavy, the projection texture coordinates are modified (displaced) along a normal vector at a given water point. This normal vector is usually obtained from a normal map. This technique can be used only for lakes, however, as these waves are way too low to represent ocean ones. An alternative to this technique is the use of dense vertices grid (the technique known as the projected grid). Vertices in this technique are transformed in the vertex shader and this way it is possible to achieve realistic waves (not an imitation as in the previous example). In the pixel shader, reflection and refraction textures are applied as previously.

A modification of this technique in turn is to use the vertex texture's fetch mechanism [7], allowing sampling of the texture in the vertex shader. Thanks to that at this stage it is possible to use a height-map. This mechanism, introduced in the 3.0 shader model, gives much better results but still they are far from being realistic and are not yet used that often for water rendering.

As you can see, the majority of existing and popular techniques focus primarily on modelling the appearance of the surface and less on the optics. I hope you remember me claiming that optics are very important.

Although these approaches are the most popular techniques for rendering water up to now (despite the ever-increasing power and the capabilities of GPUs) there are several problems associated with them:
The transition between the shore and the water is very hard and thus unrealistic. The resulting water often resembles metallic substances like mercury. In nature, water is free from sharp edges – they are always very soft and smooth. Even liquids with a very high density (like oils) do not have them.
To obtain good quality waves it is required to have an adequately dense mesh, otherwise the waves will have sharp edges.
It is difficult to obtain the correct colour extinction with depth, as well as the effect of the transition from the shallow water to the depths. This is due to the fact that usually in the process of rendering water we have no information about the depth of water at a specific point of the scene. Thus, the water is simplified, and has incorrect colour. Some implementations make use of a so-called depth map for that – a static texture describing the depth of the water area. However, even the appearance of very large objects on the bottom (for instance a ship wreck) will not change water colour whereas it should.
The solutions are rather inflexible - a share of them can be used only for rendering oceans, some for lakes.
It is recommended that some form of LOD techniques will be used in order not to process too much information about distant parts of the water from the observer.
Rendering multiple different areas within a single scene is often difficult.

Presented approach

The presented approach is based on the simple fact that water is drawn in the phase of image post-processing and is not associated with any geometry. In addition, I make an assumption of the use of deferred shading since implementation in this case is simpler and seems to be more natural. Thanks to these assumptions at this stage we have the geometry buffer filled with data. In particular, we have easy access to information about the position of each vertex visible from the camera position point of the scene. Since the geometry of the buffer is nothing more than the geometry stored in the image space, modifying it to impact the actual geometry of the scene is possible. This simple fact allows us to completely change the approach to rendering displacements and bumps. It becomes possible to move many algorithms from the geometry stage to the post-processing stage. However, I will only show how to make water this way. At the same time, we can easily get rid of the disadvantages of the traditional techniques of rendering water, which were mentioned in the previous section. The presented approach does not, however, focus on the animation of waves propagation in the water. The decision to use static height-maps, as in the example, or dynamic ones modified in accordance to the FFT, is up to you.

Modifying existing geometry

Notice that if you make modifications to the texture storing the scene point positions it becomes possible to achieve a water surface with convincing waves. We can think of this position-storing texture as of depth of water at any given point. If you know the position of the water surface L and the position P of the scene pixel, then depth = L - P. Depth = 0 corresponds to the water surface and depths < 0 correspond to the points located above the surface and so they can be skipped as they do not require further processing.

In order to obtain waves, a traditional height-map in greyscale is used to alter this depth.

The algorithm to create waves is relatively simple. It relies on tracing the ray from the position P of the scene to the observer and extruding waves in this direction. Several iterations are done. In each of them we sample our height-map and bias the water surface level by this value.

For each scene point P:
Current level of water surface L = level of water surface.
If P.y > L + H (maximum wave height), then end, because point P is above the water surface.
Calculate the eye vector E as a difference between the pixel P and the observer position: E = P – Observer Position.
Normalize E: E = Normalize (E)
For the current level of water surface L, find a point of intersection with the vector E. Mark it as the position of water point S. In other words, find the point of intersection of the plane W = (0, 1, 0,-L) with the vector E.
Perform n iterations
Sample the height-map at the point defined by S and in the direction of the vector E. The result is bias B.
Multiply B by the maximum wave height H: B = B * H
Set new L as: L = L + B
Find new S on the basis of new L value
Calculate the amount of accumulated water along the ray as A = length (P - S)
Calculate the depth of water at this point as D = S.y - P.y

Note that the results, especially the amount of accumulated water A and depth D will be used later to simulate several aspects of water optics.

The computation of normal vectors

To have realistic and convincing water it is essential to compute normal vectors. While in the case of a vertices grid, normal vectors are known already at the stage of processing geometry in the vertex shader, in the case of the presented technique, they must be calculated entirely in the pixel shader as there is no real geometry. Luckily, a simplified way of computing normal vectors known from terrain rendering applies here as well. In order to calculate normal vectors, the height-map has to be sampled in four adjacent points to the processed one, that is:

N = {W - E, 2d, S – N}

W, E, S and N are sampled directly from our height-map.

Although these normal vectors could serve as the ones used in the lighting and shading calculations, much better results can be achieved using an additional normal map. This is due to the fact that normal vectors so far are way too smooth whereas to achieve good water quality it is extremely important that detail will be present. It can be obtained only by using normal vectors with much higher resolution.

To do this we are going to use the traditional normal mapping technique. Since there is no geometry, we do not have information not only about the normal vectors (which are already set) but also about the tangent and binormal ones which are necessary for the validity of the calculations. There is an approximate method of their computation in the pixel shader, fully described by Schuler [8]. Thanks to that, the construction of the matrix necessary to carry out normal mapping is possible and pretty easy. The code is as follows:

float3x3 compute_tangent_frame(float3 Normal, float3 View, float2 UV) { float3 dp1 = ddx(View); float3 dp2 = ddy(View); float2 duv1 = ddx(UV); float2 duv2 = ddy(UV); float3x3 M = float3x3(dp1, dp2, cross(dp1, dp2)); float2x3 inverseM = float2x3(cross(M[1], M[2]), cross(M[2], M[0])); float3 Tangent = mul(float2(duv1.x, duv2.x), inverseM); float3 Binormal = mul(float2(duv1.y, duv2.y), inverseM); return float3x3(normalize(Tangent), normalize(Binormal), Normal); }

where:
Normal - normal vector, in our case, the vector set a moment ago.
Position - position in the world space,
UV - texture coordinates.
Having this matrix we can now sample our normal map texture in a traditional way. The overhead is only a dozen or so arithmetic instructions associated with matrix construction above but believe me the result is worth its cost:

float3x3 tangentFrame = compute_tangent_frame(normal, eyeVecNorm, texCoord); float3 normal = normalize(mul(2.0f * tex2D(normalMap, texCoord) - 1.0f, tangentFrame));

Just be aware that the normal vectors have to change over time. If not, the water will resemble a rigid material like plastic. Therefore, the normal map has to be sampled several times with texture coordinates varying over time. This way it is possible to achieve interfering waves of different sizes that really looks fantastic.

Optics

As already mentioned several times in this article, optics play a key role in the convincing appearance of water. In this section we are going to focus on its individual aspects.

Reflection and refraction of light

The proposed technique does not support any new way of rendering the reflections. Thus, it should be done the traditional way. This means that the whole scene must be rendered to the texture with an altered view matrix (or world one). Picture below presents the concept:

The idea of rendering the scene to the reflection texture

The location of the observer is reflected with respect to the water surface. The forward vector is flipped as well as the up vector has to be recomputed to match the others.

Water surface is also set as the user clipping plane to avoid rendering geometry above it as it would cause the reflection texture to contain too much data and therefore being invalid. For DirectX 10 and newer, SV_ClipDistance semantic should be used instead. Performing this step at the pixel shader level is not the very best idea as it is just too late – the pixel shader would be run for every pixel and thus more operations would have to be performed than is really necessary. Even if a pixel will be discarded at the very beginning of the pixel shader.

In contrast, in the case of refraction we can make some simplifications, which for most users will remain negligible. Many of the post-process effects use information about the current state of the frame buffer as they modify or operate on it. In the case of water which is treated in the paper as such we can also benefit from this. This way you will not be rendering the scene to the next texture, but even so the final result will remain satisfactory. This solution also better fits the idea of deferred shading i.e. not rendering the same geometry many times. So what has to be done is to put frame buffer content on the screen and modify it slightly to create the impression of movement with the waves. In my implementation I have just changed screen space quad texture coordinates using time and the sine function.

Specular

Another important factor affecting the quality of the water effect is specular highlighting, in some implementations also called glare or sun glow. Water is characterized by high shininess. In this article we take into account only specular caused by global light – sunlight. Local lights influence water in lesser extent and therefore they can be skipped without sacrificing too much quality. The calculation of sun glare may be done in a number of ways. In my opinion, the best results can be achieved using this snippet I found some time ago:

half3 mirrorEye = (2.0 * dot(eyeVecNorm, normal) * normal - eyeVecNorm); half dotSpec = saturate(dot(mirrorEye.xyz, -lightDir) * 0.5 + 0.5); specular = (1.0 - fresnel) * saturate(-lightDir.y) * ((pow(dotSpec, 512.0)) * (shininess * 1.8 + 0.2)); specular += specular * 25 * saturate(shininess - 0.05);

The key here is the first line of code, which reflects the eye vector so that the incidence angle is equal to the emergent angle. Therefore, an angle between normal and normalized eye vector is found. In the next few lines there is only a slightly modified process of specular calculation. The constants’ values can be changed but after testing several ones I think these gives the best results. For the shininess parameter I suggest values in the range 0.5 to 0.7.

Colour extinction

In many implementations of the water effect, light extinction is ignored. If it is implemented it is often simplified to multiplying water colour by some bluish shade. However, as I wrote in the section “Theory behind water” it is one of the most important factors affecting the apparent colour of the water. In the proposed solution light extinction is divided into two phenomena:
Colour extinction with a depth that makes all objects have bluish colour from a certain depth onwards
Colour extinction with increasing distance from the observer, the so-called horizontal transparency of water.
Depending on the depth at the given point and the distance from the observer, water will have a different colour. We first define the vector of colour extinction, which is responsible for the rate of extinction of the r, g, b components of light. On the basis of the table from the theory section this vector can look this way:

extinction = [4.5;75.0;300.0]

Then the first attempt to compute proper water colour can be as follows:

where:
refraction – colour of refracted scene,
D - the depth value determined at the end of “Modifying existing geometry” section
Unfortunately this formula has a few flaws:
The colour fades away only with depth, however there is no decay along the eye vector. This makes water pixels far from the observer have the same colour as ones near it, whereas they should be much darker.
Moreover, control over water colour is limited and very difficult. While the clean waters are relatively easy to obtain, muddy waters like lakes are much harder. It is because in the above formula, water colour is not used explicitly.
The above calculation does not take into account the true colour of water, which is bluish. In this model it is white.
We have to therefore modify our computations:

We introduced several new quantities to our formula:
The surface water colour (surfaceColor), which can be treated as a true colour of the water
The colour of deep water (depthColor)
Amount of accumulated water A as mentioned in the section “Modifying existing geometry”
Introduction of the horizontal visibility (visibility). The smaller the value of this parameter, the less transparent water will be.
In the particular case for surfaceColor = 1 and depthColor = 0 final colour of the water will be very similar to that from the previous solution. What we get is the final refraction colour. Then using the value from Fresnel term we blend reflection with refraction and add specular to the result. This way, however, water will have hard shores and that was what we strove to avoid. Therefore, we once more blend the result with refraction to an extent determined by the input parameter specifying shore hardness (1 by default) multiplied by parameter A (water accumulation). By using the shore hardness parameter we will still be able to obtain the hard edges (sometimes they can still be useful, especially when rendering NPR - non-photorealistic rendering).

Possible improvements

One of these add-ons that you can get with a small amount of work and yet which significantly improves the realism of the water effect is foam. It arises when the water hits the shore and at the tops of more choppy waves. The first of them, i.e. coastal foam, can be obtained by making an assumption that foam begins at the edge (that is, at the depth equal to 0), and to a certain depth H₁ remains constant, and in the range extincts completely.

Whereas the latter can be achieved by using the following formula:

where:
H – current height,
H₀ - height at which foam appears
H_max - height at which foam dies out.
This gives you information how much foam is visible. However, to really get it on the screen you have to sample some foam texture and multiply its colour you get by the value found above. Using a photo of coastal foam works good for that. Another possible extension is to add the interaction between world objects and water. The basic form of interaction will be provided if you use foam as described above. Since it is closely related to the depth at the point of the water area, it will appear wherever the depth is small. As a result it will appear whenever an object falls into the water.

However, the possibilities of interaction between water and the rest of the scene are literally unlimited. By modifying height-map either on CPU or in pixel shader you can easily introduce local disturbance for instance.

Even though I consider foam as an “improvement” it greatly improves general image quality and makes it warmer. Therefore you should regard using it as often as possible.


Disadvantages of the presented technique

Although the presented technique removes most of the defects of existing water rendering techniques listed in the section “Traditional approaches to water rendering” there is still enough room for improvements due to several drawbacks of the current algorithm:
It increases fill-rate as it is done in post-process and at the same time based on deferred shading. Fortunately, in most cases the water should not occupy more than half of the screen. On the other hand, often there is no need for LOD, since most of the screen pixels will be near the observer.
Water looks not that good at oblique angles. In order to get rid of this problem, it is possible to apply a different bump mapping technique. The alternative is to move the normal vector towards the observer for distant waves, which will improve their quality, as their back faces will be invisible.
In the presented approach local lights do not affect the colour of water. You can try to move water rendering before the lighting phase (we can then modify the normal data). This, however, will require modification to the existing rendering chain.

Final results and summary

You can achieve many different kinds of water using the presented technique, from tropical warm seas to cold lakes.

The algorithm presented in this paper is a good alternative to the popular techniques used for rendering water effects. By using calculations based on the observation of the phenomenon, rather than on existing mathematical models, I have managed to achieve a realistic effect, which in the case of properly selected parameters’ values makes it possible to achieve results comparable with much more expensive models. It is also a flexible solution - the appearance of water can be controlled through a number of attributes. This makes it suitable for a wide range of effects no matter it be lakes or oceans. You can even have several areas of radically different waters at the very same time.

Besides, as it does not rely on geometry it eliminates the most common flaws of the popular techniques discussed in more detail in the section “Traditional approaches to water rendering” of this article.

A word on implementation

The provided shader has been implemented using HLSL and DirectX 9. At present it requires a SM 3.0 capable GPU to run, however, it can be heavily optimized as my main concern was to make it as readable and simple as possible. Several textures are used in the code:
heightMap – height-map used for waves generation as described in the section “Modifying existing geometry”
backBufferMap – current contents of the back buffer
positionMap – texture storing scene position vectors
normalMap – texture storing normal vectors for normal mapping as described in the section “The computation of normal vectors”
foamMap – texture containing foam – in my case it is a photo of foam converted to greyscale
reflectionMap – texture containing reflections rendered as described in the section “Reflection and refraction of light”
Note that positionMap is part of the G-Buffer as deferred shading is used in my implementation. In the case of forward rendering it should be depth map. positionMap in my case stores data in view space whereas in the accompanying shader all calculations are done in world space to simplify things a bit. Therefore position data has to be multiplied by the inverse view matrix. As water is rendered as a post-process effect, a full-screen quad has to be rendered on the screen with the water shader applied.

Further reading

[1] Toman W., "Deferred shading as an effective lighting technique”, IGK’2008
[2] Placeres F. P., “Fast Per-Pixel Lightning with Many Lights”, “Game Programming Gems 6”, Charles River Media, 2006, Boston
[3] Hargreaves S., “Deferred Shading”, GDC, 2004, available on-line
[4] Guillot B., “A reappraisal of what we have learnt during three decades of computer simulations on water”
[5] Jensen L. S., Golias R., "Deep-Water Animation and Rendering“, available on-line http://www.gamasutra.com/gdce/2001/jensen/jensen_01.htm
[6] Johanson C., „Real-time water rendering - projected grid concept”, available on-line http://graphics.cs.lth.se/theses/projects/projgrid/
[7] Kryachko Y. „Using Vertex Texture Displacement for Realistic Water Rendering”, GPU Gems 2
[8] Schuler C., “Normal mapping without precomputed tangents”, ShaderX 5

Moving Noise

Wed, 27 May 2009 10:07:26 +0000
Introduction
In July 1985, Ken Perlin introduced his much-used algorithm for producing organic randomness, now known as Perlin noise [1]. That function consists of four important subroutines:
Random gradient vector lookup.
S-curve evaluation.
Grid-point dot product evaluation.
Interpolation across all dimensions.
One can find out the basics of this algorithm from the oft-visited tutorial website "Making Noise" [2]. Assuming that you are familiar with the Perlin noise function, you may have been excited to see the 2002 Ken Perlin paper "Improving Noise" [3]. That brief report on "improved noise" gave a very elegant implementation of the classic noise algorithm useful to so many procedural methods. Other optimizations and extensions have been made to the Perlin noise algorithm for use in various settings. Among these is "simplex noise" [4][5], offered by Perlin in 2001 - a straightforward approach to quickly computing high-dimensional noise functions - and "modified noise" [6], by Marc Olano - which produces certain desirable noise properties efficiently on GPU systems.
Graphics using the 2002 improved noise algorithm (D Nielsen)

  I developed the following "dynamic noise" algorithm in March 2004. In this context, the word "dynamic" refers to the inherently mutable nature of this noise function. In other words, it moves. While rereading my senior design project from the year before (a hardware-description-language implementation of improved noise with minor novelties), I noticed I had obliquely mentioned that the noise algorithm could likely be made dynamic. With a bit of tinkering, this algorithm emerged. I do not know whether something similar has arisen any place in the field, but the function has proven useful, so it seemed time to make certain it was passed on. (All other attempts to do so have seemed unusually thwarted by the universe.) Okay, enough about myself.

The general enhancement of this algorithm is noise movement without need of slicing through a higher dimension, so much less computation is required. The general drawbacks are slightly augmented memory usage as compared to traditional noise due to scaling of the gradient vectors, and use of an iterative timestepping approach that makes backtracking to a previous noise state or changing flow velocity more difficult than in the traditional approach of slicing through a higher dimension.

Method
We should begin by describing a 2D implementation of improved noise. Perlin provides code for the 3D case accompanying his 2002 brief, but leaves other cases out. We can use our own sensibilities to create a 2D case. The main difference is the set of possible gradient direction vectors. We will place them within the set {[1,2], [1,-2], [-1,2], [-1,-2], [2,1], [2,-1], [-2,1], [-2,-1]}. These vectors do not align with either the grid lines or the grid diagonals, as Perlin suggests they should not.



  Generating dynamic noise is very similar to making improved noise, with one primary difference: all 3-bit gradient direction vectors are scaled by an unsigned value. This produces a slightly less dense noise function, but a more controllable one. For the 2D case, we will scale the gradient vectors with a 3-bit unsigned value. The 3-bit [0..7] range simply emerged from guesswork and testing. It seems to be the smallest bit-width that produces mostly decent results, at least with a constant change factor in scale. The change in gradient vector scale per time step is stored in a 1-bit value, representing +1 or -1 change. When changing, the gradient vectors just scale up and down, up and down: 0, 1, 2, 3, 4, 5, 6, 7, 6, 5, 4, 3, 2, 1, 0, 1, 2, ... The gradient vector scale change is made based on a coin toss: true then change, else leave it alone. I believe I had also tried other nonlinear methods like bit-shifting scale and changing the scale sinusoidally to get higher-order continuity over time, but they did not work out so well.

All this yields the following values in the gradient LUT for the 2D case:

3 bits for scale, indicating a value in range [0..7].
3 bits for direction, indicating [1,2], [1,-2], [-1,2], [-1,-2], [2,1], [2,-1], [-2,1], or [-2,-1].
1 bit for scale delta, indicating +1 or -1.
For the 3D case, four bits are needed for direction, giving 4b (direction) + 3b (scale) + 1b (scale change) = 1B. The magic is that, when gradient vectors are scaled by zero, we can swap around the direction vectors without any noticeable effect. Then, when they are rescaled, the noise function takes on the new form determined by the gradient direction vectors. As gradient vectors change scale, new vectors are always reaching zero scale, and so gradient vector swapping occurs every time step. If you are familiar with the improved noise function, this may sound pretty simple. And it is!

This approach also implies the constraint that gradient direction vectors must be stored in writable memory, of course, not constants.

The 1D case is trivial, and the 3D case can be extrapolated from the Perlin improved noise code and an understanding of the dynamic noise approach. Higher dimensions are not presented.


Graphics using dynamic noise (D Nielsen)


Economy
The dynamic noise implementation presented here requires the following measures over static noise:

The table holding permutation values must be independent of that holding gradient lookup values.
The dot products evaluated at grid points must be scaled.
Time changes iteratively, not as a function parameter.
Gradient values must be updated whenever a time step occurs.
The true cost of these new requirements is small and is further diminished due to constraints and provisions already inherent in target systems. For example, a hardware implementation of the static noise function is well-suited to pipelining; therefore, it may be advantageous to place the permutation table in a processing stage independent of the table holding gradient vectors, in order to load-balance stages in the pipeline. Notice also that the gradient direction values require fewer bits per entry than do the permutation values, so that a separate gradient direction table does not require much additional memory.

As an approximation, we can think of the dynamic noise function as being half the computational complexity of a comparable improved noise function. A static improved noise function on a regular grid in N-space is usually computed as an interpolation of two noise functions in (N-1)-space (except with grid-point evaluations using vectors in N-space). By reducing the dimensional requirement by one dimension, less than half of the computations are required for a dynamic noise evaluation.

It is expected that dynamic noise benefits from GPU vector hardware operations, because such operations could be used to perform the grid-point dot products and scaling, and that matrix operations might be used to perform groups of vector operations for even better savings.

Conclusion
This algorithm is quite simple and useful, and can produce convincing results for most noisy effects, although the particular implementation given here is likely not suitable for very highend visuals due to coarse timestepping and efficient integer math. (For a critique of Perlin noise in general, see [7].) The believed enhancements that dynamic noise provides are reduced computational complexity by over half, better localized noise density conservation (so there is more of a sense of flow), and no question of slice orientation or ordered holes appearing from use of a higher dimension. This method also provides smoothly variable flow velocity when rotating a slice through a higher dimension. While the need for independence between the gradient table and the permutation table is required in this dynamic noise implementation (unlike the reusable table of static noise), a pipelined hardware implementation begs this independence anyway.

A set of dynamic noise C++ classes (1D, 2D, and 3D) can be found in source code below, along with a separate Windows demonstration binary. The executable displays layered dynamic noise in various places around a scene. These files are provided for your convenience. There is no guarantee about any behavior on any system - use and modification is allowed within your own liability. noise.zip

Thanks and happy coding!

References
[1] Perlin, K. An image synthesizer. 1985.
http://portal.acm.org/citation.cfm?id=325247

[2] Perlin, K. Making noise. 2000.
http://www.noisemachine.com/talk1/index.html

[3] Perlin, K. Improving noise. 2002.
http://mrl.nyu.edu/~perlin/paper445.pdf

[4] Perlin, K. A sheet of simplex noise. 2006.
http://mrl.nyu.edu/~perlin/homepage2006/simplex_noise/index.html

[5] Gustavson, S. Simplex noise demystified. 2005.
http://webstaff.itn.liu.se/~stegu/simplexnoise/simplexnoise.pdf

[6] Olano, M. Modified noise for evaluation on graphics hardware. 2005.
http://www.cs.umbc.edu/~olano/papers/mNoise.pdf

[7] Cook, RL, et al. Wavelet noise. 2005.
http://portal.acm.org/citation.cfm?id=1073204.1073264

Dynamic 3D Scene Graphs

Wed, 10 Dec 2008 10:00:08 +0000

Introduction
A scene graph is a very common structure used in computer graphics to organize the visible objects of a scene. By maintaining a hierarchy of objects in a scene, it becomes very intuitive to work with an object and its children. For example, if a vehicle object is added to a scene, whose children are its wheels, then by moving the vehicle we can automatically move the wheels as well. This automatic movement of the wheels is possible because the wheels are positioned relative to their parent, the vehicle. Relative positioning is a powerful concept, and is fundamental to the infinite scene graphs discussed in this article. Essentially, we must remember that whenever we position something in a 3D space, it is always positioned relative to something else, whether that something else is a vehicle, a planet, or a game universe.

Most scene graphs used in computer games today are restricted to being scene trees, as there will be a root scene node, and multiple generation of children underneath that, but never any cycles. For example, a common situation might be to have a city as the root node, with buildings and roads as its children. The children of the buildings will be the rooms, and the children of the rooms will be furniture. It is very uncommon in computer game scene graphs to see a scene graph node referencing itself or one of its parents, as it is more difficult to see any immediate application of how this would be useful. In our city, an example of a scene node referencing its parents is if, on the table inside a room in the city, there is a snow globe with a city inside of it, whose children are roads and buildings and etc...

Extending a scene tree to a scene graph allows us to create fractal-like objects with infinite details. In games, this means objects like trees or plants can be defined that have infinite detail. It allows the construction of unique worlds that would have otherwise been impossible to construct, such as worlds with the world contained in itself all over again. Some of the implementation details required for a true scene graph to be feasible will also be useful in their own right, allowing for massive scenes including objects of vastly differing magnitudes. Imagine zooming out from a leaf, to a tree, to a landscape, to a planet, to a solar system, and then to the entire galaxy, seamlessly.

This article will explain how to achieve an implementation of a true scene graph, as opposed to a scene tree. It will discuss how the scene graph idea itself, as well as the implementation details necessary to support it, can be applied to solve common and not so common problems in video games and other domains.

Defining the Scene Graph
In order to begin discussing implementation details, we must first determine how to describe a scene graph. Imagine a man with his arms out, holding in each hand, another man with his arms out, holding in each hand, another man with his arms out, ad infinitum. If we were to explore a 3D scene like this, it might appear that there is an infinite amount of detail, so how can we hope to describe all of it? This is not quite true though. Just like how a recursive function is defined in a finite amount of code, we can also define this man scene in a finite amount of information.


Figure 1: A scene graph where a node references itself

  First, let's formalize the information in a few nodes of a standard scene tree. Imagine a scene consisting of a garage, with a car inside of that, and four wheels on the car. We might write this as follows:

Garage -> (GarageModel, M₁), (Car, M₂)
Car -> (CarModel, M₃), (Wheel, M₄), (Wheel, M₅), (Wheel, M₆), (Wheel, M₇)
Wheel -> (WheelModel, M₈)

Where the line “X -> (Y₁, M₁), (Y₂, M₂), ...” reads: for all i, Y_i is a child of X, positioned relative to X by the linear transformation represented by the matrix M_i. In our example, Garage, Car and Wheel are the internal scene tree nodes, and GarageModel, CarModel and WheelModel are the leaf scene graph nodes, or terminal nodes. The terminal nodes have no children, and these are the nodes that actually result in a model being rendered when traversed during the scene graph pass.

In other words, drawing a garage involves drawing a GarageModel and a Car (which involves drawing a CarModel and four Wheels (which involves drawing a WheelModel)). Drawing any of the terminal nodes, in this example, implies rendering the model of the object (IE: a car mesh loaded in from 3D Studio Max).

Now, in order to allow for the definition of a true scene graph as opposed to a scene tree, we simply allow any of the Y_i's in the line “X -> (Y₁, M₁), (Y₂, M₂), ...” to be X itself, or a parent of X. We discussed the example of a man holding himself above, and now we shall formally describe it:

Man -> (ManModel, I), (Man, M₁), (Man, M₂)

Where I is the identity matrix, and M₁ and M₂ are linear transformations placing the unit cube in to the ManModel's hands. Figure 2 shows a graphical representation of the above, Man definition, as seen in the provided implementation tool.


Figure 2: A graphical view of the "Man holding himself" scene definition. In the SceneComposer sample application

  We now see that this seemingly complicated scene is actually very simply described.

Implementation Details
Now for the fun part. We will first discuss the difference between the Scene Graph definition, which is described in the section above, and the Scene Graph instance, which is a tree. We will then discuss how to terminate the recursion so that the loops introduced by a scene node referencing itself or its parent will not cause our rendering loop to go on forever. Finally, we will discuss how to maintain smooth transitions from very large objects, such as a galaxy, to very small objects, such as leaves on a tree.

Scene Definition vs. Scene Instance
Before moving on, we should discuss the difference between a scene graph definition, and a scene graph instance. While the definition of a scene graph is finite (for example, the “man holding himself” scene definition consists of only one line), the instance of a scene graph used for rendering will not be. The scene graph instance will always be a tree. The instance of a scene graph will be the tree formed by the following algorithm, whose input is the node in the scene definition to start instancing the scene at:

  SceneInstanceNode CreateSceneGraphInstance(SceneDefinitionNode instanceStart):         // Instance the current node    SceneInstanceNode instancedNode = new SceneInstanceNode(instanceStart); // Instance the children of the current node    ForEach child of instanceStart:        instancedNode.AddChild(          CreateSceneGraphInstance(GetSceneDefinitionNode(child)),          GetMatrix(child));    return instancedNode;  The application of this algorithm on the “Man holding himself” scene definition is shown graphically in figure 3. The transformation from a scene definition to a scene instance is important, because since the scene instance is a tree, we can now simply process (IE: render) it using standard “scene graph” processing techniques. Obviously, if the Scene Definition graph has loops in it, then the algorithm presented here will produce an infinite scene instance. The title of this article gets its name from the fact that the scene tree which results from the instancing process is dynamically constructed.
Figure 3: Example of CreateSceneGraphInstance algorithm applied to "Man holding man" scene

  Instancing of the scene definition can be implicit or explicit. Explicit scene instantiation is described above, and involves formally processing the scene definition in to a scene instance, and then processing the resulting scene instance like a normal scene tree. Implicit scene instantiation would be to process the scene definition as it is traversed. That is, no “Scene Instance” data is explicitly created.

Explicit scene instancing can allow for a more powerful and robust system, since every instantiated scene node can now be treated as its own object. That is, if a person scene node happens to be instanced twice, each of these instances can now maintain their own, separate properties, such as animation and state, even though they originated from only one scene definition node. The downside to explicit instantiation though, is that it requires more memory.

Terminating the Scene Graph Instancing
As with all recursive algorithms, there must be a base case or else we will recurse forever. By restricting the linear transformations of the children of a node in the scene graph to be contractions (that is, the transformation has a component that scales the object down in size), we ensure that all of the children of a given node are smaller than that node. This implies that if we traverse down the nodes of a scene graph, we will arrive at smaller and smaller nodes. So, to terminate our scene instancing algorithm, we should stop recursing when our scene nodes are “small enough”.

So when is a scene node “small enough” to stop recursing? A good method for determining this is to use the screen-space size of the node about to be rendered. This can easily be approximated by determining the maximum screen-space size of the bounding box for the given scene node. In order for us to calculate the screen-space size, we must have access to the model-perspective transform, and so this must be passed down our scene instancing algorithm. At each node, we should apply this transformation to that node's bounding-box, giving the screen-space points of the bounding-box. Since an axis-aligned bounding box is easier to deal with, we create a 2D axis aligned bounding rectangle enclosing the 8 screen-space points of the cube. Now it is simply a matter of determining the area of this axis aligned bounding rectangle (IE: width * height) to give an upper bound on the screen-space area of the node.

With the screen-space area of a scene node available to us, we can check whether it is below a rendering threshold or not. If the screen-space area of the scene node is too small, or below the threshold, then we return early, and this particular node, or any of its (possibly infinite) children do not end up in the scene instance tree. Since each recursion of the scene instancing algorithm takes us to smaller and smaller nodes, we are guaranteed to eventually arrive at nodes that are below the instancing threshold, guaranteeing that the algorithm will always terminate. What we are left with is only a part of the (possibly infinite) scene instance tree, but this is the part that is most relevant to us, since everything that was cut are scene nodes that occupy very little screen space.

Frustum culling can also be effectively applied at this stage to farther reduce the size of the resulting instance tree. This is easily done since we already calculate the screen-space bounding cube points of each node.

Infinite Zoom

The Zooming Problems
While enough has been explained already to produce a working true scene graph rendering system, it's usefulness is currently limited to allowing for objects whose level of detail is increased the closer you come to them (for example, a plant defined recursively will come in to more detail the more screen-space it occupies). While this is already a nice improvement to existing scene trees, we can make things a lot more interesting by allowing for the user to zoom infinitely closer to items in the scene.

Our first attempt might be to allow zooming by simply decreasing the camera's field of view. Since we are passing down the model-projection matrix to our scene instancing algorithm, the lowered field of view will be taken in to account and the level of detail will be adjusted accordingly. Unfortunately, it becomes very difficult to navigate around an object of focus we have zoomed in to because the slightest rotation of the camera will target a completely different space. Also, we run into floating point accuracy problems when zooming in too much.

Stepping back from the problem a bit (no pun intended), we realize that another way to achieve “zooming” is to simply move the camera closer to the object. Due to the mechanics of the the perspective projection matrix, this has the effect of making objects bigger (or smaller when you move away from them). Using this method allows us to maintain the same field of view, so camera rotations have the same degree of effect when “zoomed in”. This method still leaves us with two problems though. The first problem is that if we continue to move our camera with the same speed we always had, it will now be moving much faster, relative to our new, smaller, object of focus. This is a problem because moving our camera 1 unit relative to a galaxy will move our camera an enormous number of units relative to a leaf, so clearly we must find a way to manipulate the speed of our camera. The second problem is that we still suffer from floating point accuracy problems, since floating point numbers do not allow us to represent an infinite amount of zoom to an object.

The Zooming Solution
The problem of floating point accuracy, and the problem of unsuitable camera speed when viewing smaller objects have a common solution. We simply need to realize that we can render the exact same scene under completely different coordinate systems, as long as everything is still positioned relative to each other. The observer of your world won't know if a leaf is 1 unit big or 100 units big, as long as it is still fifty times as small as the tree it is attached to. What we need to do is be able to dynamically select different scene nodes as being the scene's reference node, the node that every other node, as well as the current model-view transformation, is relative to.

An Alice in Wonderland game, for example, might start with the world being rendered with the frame of reference being a large scene containing a landscape with a rabbit hole in it, and all of the geometry consisting of the interior of the rabbit hole, with a desk at the bottom of it that has a little bottle with the words “DRINK ME” on it. Within this scene graph might also be a much tinier room with a tiny, locked door in front of it. This smaller room will have its own scene graph children in it as well. When Alice drinks the contents of the bottle, the game engine can then make the smaller room the new reference, everything else being made relative to it. As long as the new frame of reference is not scaled on an entirely different order of magnitude from the old frame of reference, the viewer will not notice when the reference frame switch is made, as the camera's position and orientation will also have changed to be relative to the new frame of reference.

By taking the transformation that positions an object relative to the current frame of reference (calculated by accumulating the matrices in the scene graph leading from the frame of reference to the object in question), and inverting it, we obtain a transformation that positions the current frame of reference relative to the object in question. The following pseudocode shows how to calculate the matrix used to change the frame of reference:

  Matrix FindRelationMatrix(sourceNode, destNode):    // Helper function for getting the list of all a node's ancestors    // In the returned list, the original node will be the first item    // and the root will be the last    NodeList GetAncestorsList(node):   ancestorList = [];   curNode = node;   while curNode:         ancestorList.Append(curNode);         curNode = curNode.parent;    // Accumulate the tranform matrices leading from the source node to    // the destination node    Matrix sourceToRoot = Identity;    Matrix rootToDest = Identity;    For node in GetAncestorList(sourceNode):   sourceToRoot = sourceToRoot * node.GetTransform().Inverse();    For node in GetAncestorList(destNode):   rootToDest = node.GetTransform() * rootToDest;    return rootToDest * sourceToRoot;  By applying this transformation to the current model-view transformation, we obtain a new view transformation that is relative to the object we wish to re-focus on. If the new frame of reference is smaller than the old one, we have just made the entire universe bigger, but since the view matrix has been adjusted to account for this as well, the universe appears unchanged. Now, though, we can begin to show even smaller objects that would have been problematic due to floating point accuracy in the previous frame of reference. By continually changing the frame of reference, we can continue to zoom in (or out) of objects infinitely. Notice also that since the view matrix is continually being modified, then so is the velocity of the camera, so that when we switch to a smaller frame of reference, the camera automatically begins to move slower. This is desirable, since it means the camera will now move at a reasonable speed to allow the viewer to analyze the detail of a much smaller object.

Rendering Scenes With Big, Distant, Detailed Objects
When viewing a very small object up close, it is possible for the viewer to be also looking at a very large object in the background. Think of the camera looking at the tip of a leaf, while behind the leaf is a planet, and behind the planet is a sun. Even though the planet and sun are far away, they're also really big, so they still occupy considerable screen-space area. We can't render a leaf, a planet, and a sun in the same scene though, because the z-buffer would provide very poor quality for the distant objects, resulting in ugly z-fighting artifacts.

A solution to this problem exists by rendering the the scene instance tree in a breadth first fashion, starting at the root. This has the effect of drawing the largest objects first. A “rendering reference frame” is created, starting with the root (biggest) node. The camera is re-positioned relative to this reference frame, and rendering begins. When recursing down towards the current object of reference (IE: towards the viewer), we track how much the scale has changed from the root node. When it is detected that the scale has been reduced past a given threshold, we clear the current z-buffer, reset the “rendering reference frame” to the current scene node, adjust the camera according to the new reference frame, and continue rendering in to this new “layer”.

The layering approach works because we are traversing the instance scene tree towards the camera, so every time we enter a child that also contains the camera, we are rendering objects closer to the viewer. A downfall to this solution is that there will be rendering order artifacts if the bounding boxes of two scene node children overlap. In this case, it would be possible to flag that the particular node is “required to be rendered with the same z-buffer as its children”, so the z-buffer will never be cleared when jumping from that node to its children.


Figure 4: A man on in a city on a planet, with another solar system in his galaxy, all rendered in the same scene


Applications
To put what's been previously described in to some perspective, a description of how the system might be used in a game like Spore will be described. I do not know how Spore implemented their transitions from planet-side views to galactic views, I only discuss here how it might be implemented using the dynamic scene graph system.

In Spore, the player can obtain a spaceship with which they can zoom out to a view of the entire galaxy, flying around from one to the other. The player can fly to a solar system in a galaxy, and then a planet within the solar system. When the player flies towards a planet, the camera zooms in to the planet surface where the player can manage the local situation. This cannot be implemented relative to a single coordinate space, because the size of a planet surface is microscopic compared to the size of the galaxy, so there would be large floating point errors, resulting in inaccurate graphics and collisions on the planet surface.

We do not wish to generate the entire galaxy and every planet within it, since if we allow for an arbitrary amount of detail on all scales, we will run out of memory. Furthermore, in the context of a game with a changing game-state, it is less trivial to apply the dynamic scene graph generation strategy, because we want the game to never forget that we built a city on planet A, even if we're on the other side of the galaxy, and planet A is now too small to be in our scene graph instance. With these challenges overcome, we would have a reasonable framework to implement a Spore-like game within. We will discuss these topics now.

Galaxies of Arbitrary Size and Detail
In order to support galaxies of arbitrary size and detail, we need a method of only fetching the detail currently relevant to us, while ignoring the rest. Suppose the planets in the galaxy and everything on them are created procedurally. When we go to instance the solar system scene graph nodes, they should be defined in such a way that the instancer is instructed to randomly position children planet nodes. The planet nodes will then be instanced as required, where they should instruct the instancer to randomly generate a landscape and populate it with cities. The cities may themselves be child nodes which are instructed to procedurally generate buildings and streets, and so on. Since the creation processes only happened during instancing, and we only instance scene graph nodes that are visible to us, then we are able to construct only the part of the universe that we are currently concerned with.

This same technique can be extended to support artist-defined universes. In this case we can load in the art assets for the player's starting point and then, during scene graph instancing, for every node that contains an art asset, we determine its children and whether or not to instance them. If we decide it's relevant to instance an art asset's children, then we load those child assets off the disk and store them in memory. Art assets should be removed in a Least Recently Used fashion, so that memory for art that is no longer relevant is continually freed. Once again, this method can be used to achieve an artist-defined universe with arbitrary detail that is completely connected, yet automatically only the portions of that universe relevant to the player will be loaded in to memory. It should be noted here that in this case, the art assets are actually nodes of the scene definition, so the actual scene definition will also be dynamically extended as new assets are loaded. This technique may also be useful for applications like Google Earth.

In order to for these more sophisticated instancing methods to function properly, we cannot blindly re-instance all our objects from the scene definition. We must remember the objects we have instanced, and only re-instance objects that have just become relevant. This is discussed next.

Persisting Instanced Nodes
Before we proceed, it should be recognized that in a game like Spore, it would be convenient to deal with a planet and everything in it as a whole, as it is difficult to represent a spherical planet with mountains on it as a scene graph. What we can do instead is that when the planet is instanced, we can also instance the randomly generated terrain using standard terrain generation techniques as well as initialize a spatial partitioning structure to deal with cities and the critters moving around the planet. Everything instanced as part of the planet will operate within the planet's frame of reference, which is given meaning by the planet node's position within the scene graph instance.

When a planet is instanced, then, we would want to initialize many items inside of it. Furthermore, since the player can affect the state of an instanced planet, we would like to remember the changes they made to the planet's state. To accomplish this, we should mark the planet instance as “unforgettable” to the scene graph system when the player makes a change. In doing so, the planet should never be de-instanced, even if it is considered no longer relevant (at a given moment). In remembering a planet, we are forced to remember all of the planet's parents as well, otherwise the system has no way of recognizing where the planet fits in with the rest of the scene graph. This is not unreasonable on memory though, since the height of the instanced scene tree is logarithmic with respect to the amount of nodes in the tree.

It is true that a player could then create a “worst case” by jumping from planet to planet making slight changes, forcing our system to remember more changes than we have room for. This is an inevitable limitation however, as the computer is asked to literally remember more things than it has a memory for. A possible solution would be to assume that if the computer can't remember all the player's changes, then the player probably can't either, and we can apply a Least Recently Used scheme to forget changes that the user has made to the game world.

In order for instance remembering to work, our scene instancing algorithm should be augmented with logic to check that a scene element has already been instanced and if so, it should re-use the existing instance instead of re-instancing. Pseudocode is given for this style of instancing, where the algorithm will take the root of an existing scene instance tree and instance any nodes that have not yet been instanced, while deleting nodes that no longer need to exist.

  CreateSceneGraphInstance(SceneInstanceNode instanceStart):          // Check if this node's children have already been instanced    if instanceStart.HasChildren():        if IsRelevant(instanceStart):          // Simply recurse in to the children, in this case          ForEach child of instanceStart:         CreateSceneGraphInstance(child);   else:          // Children are no longer needed, get rid of them          instanceStart.DeleteChildren();    else:   // Children have not been instanced.   If IsRelevant(instanceStart):          // And they are relevant, so instance them          instanceStart.InstanceChildren();  A sample implementation is provided. It includes two programs, Scene Composer and Scene Explorer. Scene Explorer can be used to quickly view a scene created by Scene Composer. Scene Composer allows you to import Wavefront Obj model files and then compose them together to create a scene graph definition. The sample scene composer allows you to specify more than one production with the same symbol on the left side (ie: “Planet -> RedPlanet” and “Planet -> Green Planet”). In this case, the scene instancer will randomly choose one of these productions when it instances a Planet symbol. Since it is randomly choosing a production to use here, we must remember the choice so that when we instance the next frame, we do not randomly choose another production. This is implemented by keeping the instance tree available for the next frame, at which time the instancer will traverse the old instance tree, instancing new symbols only when they have not already been instanced. Symbols that become irrelevant are destroyed.

Some example scenes are provided. The “Fern” example attempts to re-create the popular 2D fern fractal, in 3D. It uses a green cylinder as the only model, representing the fern's “stem”. It then uses the fern symbol recursively to represent the fern's branches. The final scene resembles a plant with infinite detail which can be zoomed in to or out of. This example demonstrates how the infinite scene graph system can be used to create objects of infinite detail.

The “Odyssey” example features a galaxy of solar systems of planets of cities of buildings with men in them. The buildings in the city have random heights, accomplished by having the “building” production randomly choose between four building heights, when instanced. The men in the scene recursively feature galaxies in their eyes, so the viewer can zoom in or out as much as they want. This example shows how one can implement an explorable Spore-like universe (although the instanced scene nodes maintain very little state separate from the scene definition).

The sample application is packaged as a Windows installer, which includes the above examples and others. The installer will place a README file in your start menu explaining how to use the SceneExplorer and the SceneComposer.

The source code is also available as a separate download as a zip file. It is built using SCons and requires some Boost libraries. I have not tested it on any compiler besides MSVC 2005 and MSVC 2008.

Conclusion
The full system described here consists of the combination of a number of techniques which compliment each other. Each of these techniques provides for interesting applications in their own right. These techniques can be applied together, or by themselves to solve common problems that arise in video games, as well as enabling new creative game-play ideas.

With a scene graph that can contain loops, we allow objects to have an infinite amount of detail, where the amount of visible detail is determined by the screen-space area that the object occupies. Objects such as plants can be defined recursively, similar to how the famous fern fractal is defined. Applying this technique can give a plant that occupies very little detail (and therefore graphics processing time) when far away, but is continuously refined as the viewer comes closer. Applying this technique allows an artist to describe an object in as much detail is they wish, letting the run-time engine determine how much of that detail to display at any time. A fern-like object constructed with the system is shown in figure 5.


Figure 5: A view of a fern object. A green cylinder is the only model used

  Using the system that allows for a change in frame of reference, we enable scene graphs where the very large objects, such as galaxies, can co-exist in the same data structure with the very small objects, such as leaves. This is a very powerful system, since a viewer can seamlessly move from a close up of the tip of a leaf, back to a wide-shot view of an entire galaxy. As long as the frame of reference is continually switched to a bigger (or smaller) one as the user zooms out (or in), one gigantic scene can describe an entire galaxy and everything inside it with intricate detail.

These techniques can also be applied to create previously unimaginable artistic continuous environments where the user can continuously zoom in to objects as far as they desire, cycling the user back to previous environments as needed. It opens the door to creative new ideas, and previously unimaginable scenes.

References
The model of the man was created in Blender by following the tutorial from this site

Computation of Bounding Primitives on the GPU

Wed, 12 Nov 2008 09:00:39 +0000

Introduction
Finding the bounds for a set of triangles is one of the most frequently-performed computations in games. It is fundamental to collision detection, scene culling, and ray tracing. Performing it on the GPU offloads work from the CPU, and can be an order of magnitude faster. The axis-aligned bounding box, or AABB, can be represented as 6 planes:

  struct AABB {    float Left;    float Right;    float Top;    float Bottom;    float Near;    float Far; };  Perhaps a more common alternative is simply storing two points for the corners:  struct AABB {    Vector3 MinCorner;    Vector3 MaxCorner; };  CPU pseudocode to compute the bounding box over a set of points is trivial:  AABB aabb; aabb.MinCorner = (MaxFloat, MaxFloat, MaxFloat); aabb.MaxCorner = (-MaxFloat, -MaxFloat, -MaxFloat); for each vertex in buffer {    aabb.MinCorner.x = Min(vertex.x, aabb.MinCorner.x);    aabb.MinCorner.y = Min(vertex.y, aabb.MinCorner.y);    aabb.MinCorner.z = Min(vertex.z, aabb.MinCorner.z);    aabb.MaxCorner.x = Max(vertex.x, aabb.MaxCorner.x);    aabb.MaxCorner.y = Max(vertex.y, aabb.MaxCorner.y);    aabb.MaxCorner.z = Max(vertex.z, aabb.MaxCorner.z); }  If you wish to leverage D3DX, simply use the D3DXComputeBoundingBox function to perform the above operation.
Naïve GPU Approach: Leverage the Depth Buffer
The aforementioned loop is embarrassingly parallelizable and therefore well-suited to the GPU. Using a general-purpose computation language like NVIDIA’s CUDA technology would make this easy, but it is often desirable to stay within the confines of the graphics API itself. How could we use the graphics API to perform a computation like this? The divide-and-conquer philosophy tells us to approach the problem by solving its sub-problems first. It can be helpful to think about only one dimension at a time. If we can simply compute the near plane, then we should be able to compute the other five planes via transformations.

Finding the “nearest” point is something fundamental to real-time graphics. Anyone familiar with OpenGL or D3D knows that where there is 3D rendering, a depth buffer lurks beneath:

Note that we are interested in only one depth value, while the depth buffer represents many depth values. The trick is to render to a 1x1 viewport. It’s difficult to visualize rendering an entire model into a 1x1 area, so it may help to think of it as the end of a continuously shrinking viewport:


Aside: Depth Buffer Visualization with Direct3D 10
To produce a grayscale visualization of the depth buffer similar to the pictures in this article, one approach using D3D10 is the following:
Create a texture that has a color format equivalent to the depth buffer format. For example, if the depth buffer is DXGI_FORMAT_D32_FLOAT, then the texture should have a format of DXGI_FORMAT_R32_FLOAT. The texture should have the same width and height as the render target.
After rendering the scene as you normally would, un-bind the depth buffer by passing null into the last argument of ID3D10Device::OMSetRenderTargets().
Use ID3D10Device::CopyResource() to copy the values from the real depth buffer into the grayscale texture from step 1.
Bind the grayscale texture to a pixel shader resource and draw it to the screen using a full-screen quad.
An example of this procedure can be found in the sample code attached to this article.
    Since the nearest point always trumps the points that render behind it, the depth buffer of our 1x1 render target should hold the nearest Z value. Computing values for the other five planes can be achieved by simply rotating the model appropriately:

There are some caveats to using the depth buffer in this way:

Be sure to use a high-precision depth buffer. In the case of D3D10, use DXGI_FORMAT_D32_FLOAT when creating your depth texture. Remember that the viewport is only 1x1, so the performance implication of this is relatively unimportant.
Because we are rendering to a 1x1 viewport, each triangle takes up a screen area that is only a fraction of a pixel. This can lead to rendering only one interior point per triangle, which results in an incorrect bounding box. To bypass this issue, use point primitives rather than triangles. In the case of D3D10, set your primitive topology to D3D10_PRIMITIVE_TOPOLOGY_POINTLIST.
It seems cumbersome to re-render the object in 6 different orientations. One idea is to use only three rotations, then set the graphics API to use reverse testing for each “Max” plane, as depicted here:

To achieve this with D3D10, set your DepthFunc to GREATER_EQUAL rather than the default LESS. Also be sure to change your clear value appropriately. So, your effect file would change as follows:

  DepthStencilState MyDepthState {    DepthEnable = TRUE;    DepthFunc = ~~LESS~~ GREATER_EQUAL;    DepthWriteMask = ALL; };  And your clear call would change like this:  device->ClearDepthStencilView(depthView, D3D10_CLEAR_DEPTH, 1 0, 0);  Unfortunately, this approach is not much of an improvement as it still requires a total of 6 rendering passes. How can we reduce the number of passes down to 1 or 2?
Exploit the Blending Hardware
When rendering, we typically render to a color buffer, not just a depth buffer alone. Each pixel in our 1x1 render target has slots for red, green, and blue that we’ve been ignoring. The secret sauce lies in the blending hardware. Typically, blending is used to compute a simple weighted sum, as in:  αc_src+(1-α) c_dest

  However, graphics APIs typically allow more than just addition, as seen in the following table.   DirectX 9 DirectX 10 OpenGL   D3DBLENDOP_ADD D3D10_BLEND_OP_ADD GL_FUNC_ADD   D3DBLENDOP_SUBTRACT D3D10_BLEND_OP_SUBTRACT GL_FUNC_SUBTRACT   D3DBLENDOP_REVSUBTRACT D3D10_BLEND_OP_REV_SUBTRACT GL_FUNC_REVERSE_SUBTRACT   D3DBLENDOP_MIN D3D10_BLEND_OP_MIN GL_MIN   D3DBLENDOP_MAX D3D10_BLEND_OP_MAX GL_MAX   Those last two rows should catch your eye!

Before we can use these blending operations, we should make sure that we’re rendering into a floating-point render target. In D3D10, this can be done with DXGI_FORMAT_R32G32B32A32_FLOAT. Despite the fact that we’re exploiting the hardware’s blending stage, we actually will not be using the alpha channel. We specify alpha in our format only because it is typically required by the API when specifying a render target. So, by simply writing X Y Z into our R G B channels while using MIN blending for the first pass, and MAX blending for the second pass, we can compute a complete bounding box! The D3D10 HLSL for this technique follows.

  struct VS_OUTPUT {    float4 Position : SV_POSITION;    float4 Pretransformed : POSITIONT; }; VS_OUTPUT FindBoundsVS(float4 Position : POSITION) {    VS_OUTPUT Output;    Output.Position = float4(0, 0, 0, 1);    Output.Pretransformed = Position;    return Output; } float4 FindBoundsPS(VS_OUTPUT In)  : SV_TARGET { return In.Pretransformed; }  Since we’re rendering into a 1x1 viewport, the vertex shader sets the outgoing position to (0, 0, 0), then writes the pretransformed position out to the pixel shader, which simply writes it out into the floating-point render target. To read the values of our bounding box out from the render target, we’ll need to copy the values into a “staging area” that can be accessed from the CPU. Because of this copy, which is slow, it is best to actually render to a 2x1 texture, and shift the 1x1 viewport between passes, allowing us to copy only once instead of twice. Summarizing the technique for D3D10:
Create two floating-point 2x1 textures with the following parameters:

BindFlags = RENDER_TARGET, CPUAccessFlags = 0, Usage = DEFAULT
BindFlags = 0, CPUAccessFlags = READ, Usage = STAGING

Bind the output merger stage to the 2x1 render target texture.
Set the viewport to 1x1 with TopLeft at (0, 0).
Set the blending state to use BlendOp = MIN, SrcBlend = SRC_COLOR, DestBlend = DEST_COLOR.
Render the geometry using the pass-through shaders described above.
Move the viewport’s TopLeft to (1,0)
Change the BlendOp to MAX.
Render the geometry again.
Unbind the output merger stage from the 2x1 render target.
Call ID3D10Device::CopyResource() to copy data from the render target texture to the staging texture.
Call ID3D10Texture2D::Map() on the staging texture to read the values of the two corners from the color buffer.
Conclusion
If you’ve got your API state set up correctly, computing a bounding-box can be as simple as re-rendering the object into a 1x1 viewport and reading back the color values. The GPU represents an enormous amount of computing power that is now typical in every home PC. By thinking outside the box (no pun intended), it’s easy to find new and interesting ways to offload work from the CPU.

Settings	FPS	SSAO time (ms)
High (32 samples front/back)	150	3.3
Medium (16 samples front)	290	0.27
Low (8 samples front)	310	0.08