Tom Hammersley's Graphics and Games Development Blog

Polygonal Area Light Sources on DX11-Class Hardware

2013-06-05T13:00:00.001+01:00

The surge in interest in Physically-Based Rendering has led to a surge in interest in area light sources. There is, however, something of a dearth of information available for the aspiring implementer.

Introduction

The objective of polygonal area light sources is to define an arbitrary polygon in space that may emit light. I'm going to make two simplifications here: firstly, that the polygon emits a constant colour across its surface and secondly, that we are not considering shadowing.

The earliest source I was able to find on this subject is James Arvo's 1994 paper on polyhedral light sources. John Snyder also presents another excellent reference. The approach taken here is to compute a vector field such that at any point in space r, we may obtain the direction of the light. This vector field is a function of the position r, but also the polygon's vertices. This direction can then be substituted into your brdf equation for a particular surface.

Arvo's algorithm is simple: for each edge of the polygon he computes a normal named gamma, and a weighting value named theta. The normal is of the triangle formed by the polygon edge and the point r. The weight theta is the angle formed by the pair of vectors from the point r to each edge vertex. He forms a weighted sum of these normals, divides it by 2*PI (to normalise the weights returned by acos) and scales it by the light's intensity. This is all you need to perform basic area lighting.

There are a couple of twists. Firstly, you typically need to negate the computed vector to get it into a form ready to dot against your normal. Secondly, you must clip the polygon against the plane of the point being shaded. It is possible to have a valid light source that straddles the plane of the point being shaded.

Basic HLSL code

Here is the HLSL code to perform the basic calculation to obtain the light vector for a polygon at a given point. Here, I process quads. This is no algorithm restriction; this is purely to minimise buffer traffic.

// Derive quad vertices in world-space

float3 quadCentre = g_quadLights[lightIndex].pos;

float3 quadAxisX = g_quadLights[lightIndex].axisX * 0.5f;

float3 quadAxisY = g_quadLights[lightIndex].axisY * 0.5f;

float4 quadVertices[4] =

{

float4(quadCentre - quadAxisX - quadAxisY, 1.0f),

float4(quadCentre + quadAxisX - quadAxisY, 1.0f),

float4(quadCentre + quadAxisX + quadAxisY, 1.0f),

float4(quadCentre - quadAxisX + quadAxisY, 1.0f),

};

float3 lightField = 0.0f;

uint polyVertex1 = 3;

for(uint polyVertex0 = 0; polyVertex0 < 4; ++polyVertex0)

{

float3 delta0 = quadVertices[polyVertex0].xyz - worldSpacePos;

float3 delta1 = clippedPolygon[polyVertex1].xyz - worldSpacePos;

float3 edge0 = normalize(delta0);

float3 edge1 = normalize(delta1);

float theta = acos(dot(edge0, edge1));

float3 gamma = normalize(cross(delta1, delta0)); // Note you can't just cross edge0 and edge1 to save the normalize - that gets you the wrong result!

lightField += theta * gamma;

polyVertex1 = polyVertex0;

}

// Note that we don't divide by 2*PI as in the paper. That is just there to account for the sum of the theta values

// We don't need it since we're doing a normalise - it's just a uniform scale

// Summing normalised vectors doesn't necessarily give a normalised result. Also flip the vector ready for light calculations

lightField = -normalize(lightField);

Clipping

We must now deal with the problem of clipping. This is a particularly awkward operation, especially when it may need to be done per-pixel. Firstly, we must establish whether we actually need to clip. Since the quads are formed using a quad centre and two axes, this may be handled using something very similar to an OBB vs plane test:

float criticalDistance = abs(dot(worldSpaceNormal, quadAxisX)) + abs(dot(worldSpaceNormal, quadAxisY));

float planeToQuadCentreDistance = dot(worldSpacePlane, float4(quadCentre, 1.0f));

if(planeToQuadCentreDistance < -criticalDistance)

{

// Invisible

}

if(planeToQuadCentreDistance < -criticalDistance)

{

// Fully visible - unclipped

}

We then use standard Sutherland-Hodgman clipping to trim the polygon against the plane:

float4 clippedPolygon[5]; // Clipping against 1 plane can add one extra vertex to the quad

uint numClippedVertices = 0;

// Clip the polygon against the plane of the point being shaded

float vertexPlaneSide[4] =

{

sign(dot(quadVertices[0], worldSpacePlane)),

sign(dot(quadVertices[1], worldSpacePlane)),

sign(dot(quadVertices[2], worldSpacePlane)),

sign(dot(quadVertices[3], worldSpacePlane)),

};

numClippedVertices = 0;

uint previous = 3;

for(uint current = 0; current < 4; ++current)

{

float currentSide = vertexPlaneSide[current];

if(currentSide < 0.0f)

{

if(vertexPlaneSide[previous] > 0.0f)

{

float4 v1 = quadVertices[previous];

float4 v2 = quadVertices[current];

float3 delta = v2.xyz - v1.xyz;

float t = dot(worldSpacePlane, v2) / dot(worldSpacePlane.xyz, delta);

clippedPolygon[numClippedVertices++] = float4(v2.xyz - delta * t, 1.0f);

}

else

{

float4 v1 = quadVertices[current];

if(currentSide > 0.0f && vertexPlaneSide[previous] < 0.0f)

{

float4 v2 = quadVertices[previous];

float3 delta = v2.xyz - v1.xyz;

float t = dot(worldSpacePlane, v2) / dot(worldSpacePlane.xyz, delta);

clippedPolygon[numClippedVertices++] = float4(v2.xyz - delta * t, 1.0f);

}

clippedPolygon[numClippedVertices++] = v1;

}

previous = current;

}

The resulting polygon may then be used in the original light field calculation code.

Image Histogram Construction In DX11

2013-05-17T22:07:00.001+01:00

Preamble

Recently, I have been researching some next-gen rendering techniques. There's some great presentations out there by the usual suspects - Epic, Crytek, DICE, Guerilla and so on. Many of the vendors also have good sample code to explore. There's a lot of buzz about certain new techniques... but there's sometimes a little bit of a gap between the academic literature, the presentations and sample code. Sometimes its a little hard for the time-pressed practioner to google an example and jump straight in.

So, I thought I'd show a little example code to help people get started.

Introduction

Tone mapping using a histogram is becoming a trend in next-generation graphics. People are starting to experiment with moving away from the DX9/360/PS3 style "compute average luminance, apply exposure, apply film response curve" algorithms towards trying something different.

Histogram based methods aren't new. An early, commonly cited example is Ward et al '97. More recent examples are things like Duan et al 2010. The basic idea is that you take an input high dynamic range image, compute a histogram measuring the frequency of different luminances in that scene, and use that to build a response curve. The response curve maps from the scene's luminance range to your output display lumiance range.

So, in one fell swoop, you can take your HDR image and convert it to an LDR image. No repeated downsampling and cache thrashing, no incorrect half-pixel or half-texel offsets, no frame-by-frame adaptation. A simple algorithm that "just works".

This algorithm fits very naturally with compute shaders in DX11. However, not everyone out there is a DX11 expert; there are many practictioners busy shipping titles on DX9 class hardware who have not had the opportunity to experiment with DX11 yet. I shall detail the implementation of a very basic histogram based tone mapping operator that may serve as the basis for future expansion.

Running a Compute Shader Over an Image

Compute Shader Introduction

Compute shaders off a great deal of general, flexible power on modern DX11 GPUs. They are very flexible and offer immense general computation power. They are not, however, a direct replacement for C++ code; they have a specific programming model that you have to follow. Whilst you can ignore some details of this model, you will not achieve maximum performance. This programming model therefore largely dictates how you must structure your algorithm.

An individual compute shader contains the inner most operation of some loop. However, the compute shader is not executed as some kind of iterative loop. Instead, many separate invocations of the compute shader code are scheduled, dispatched and retired. Each invocation of the compute shader is known as a thread, and each thread receives unique addressing information. And these threads are executed in parallel. Threads are grouped into thread groups. Each thread is therefore separately executing some part of the desired operation. It is very akin to a parallel map operation in a functional language.

A compute shader is therefore unlike a general CPU process or job. It is unlike the jobs or task programming model of the Cell processor. It is very similar to shading a pixel. In contrast to pixel shading, compute shader threads can communicate with other concurrent threads. A thread may communicate with other threads in the same group using some on-chip fast memory known as group shared memory.

The conclusion of this is that when processing an image with a compute shader, you do not sit in a big x * y loop, fetching pixels and processing them. You do not DMA in parts of an image, process them, and DMA them out again. You write code that processes each pixel independently and dispatch that work to be processed as a whole in parallel. This requires a certain structure to your code, especially if you want to obtain maximum parallelism.

Tiling

A common paradigm when executing a compute shader over an image is to break the image up into tiles. Each thread typically corresponds to an input pixel.A thread group maps to a tile. Each thread is executed and aggregates data to the local, on-chip memory known as group shared memory. When all the threads in a group have completed, the results are exported from the transient group shared memory to some persistent off-chip memory.

But what is the reasoning behind this model?

If we liked, we could dispatch a compute shader with a thread group size of 1 x 1 x 1 thread. Inside this one thread of one shader, we would iterate through every single pixel in an image and perform our computation. But, this model offers very little parallelism. We would only enjoy the SIMD instruction level parallelism, and would have no thread-level parallelism. GPUs are built to be very, very wide. We would only be occupying a fraction of the GPU.

Ideally. we want to launch as many threads as possible. This undermines any kind of approach where one thread iterates across many pixels. We want a thread per pixel.

And this is what gives rise to the image tiling model. Here, we break the image into tiles of 16x16, 32x32 or whatever size yields the maximum threads and therefore parallelism possible. We exploit group shared memory rather than work directly in off-chip memory for performance reasons. A modern GPU will not simply sit and churn through every every thread of your compute shader until completion. It will execute wavefronts (groups of threads) from multiple different shaders simultaneously across the GPU. If we can read from off-chip memory but write to on-chip memory, we minimise our bandwidth consumption and play better with other shaders. Ultimately we will need to export the results, but we can do that once in a burst rather than constantly causing memory hazards.

We may need to execute subsequent passes that combine together this intermediate metadata to yield a final result. This arrangement is akin to the common map-reduce parallel programming model.

I will now dive straight into the HLSL compute shader code. For clarity and brevity, I won't present the detailed DX11 CPU-side code; I shall just detail what it needs to do.

Computing the Per-Tile Histogram.

Our first task is to compute a histogram for each tile of the image. These will later be combined to form a histogram for the whole image. Each histogram will map the scene's luminance range across N bins. Each bin holds a uint which represents the number of pixels in that luminance band.

To perform this, we need a few things on the DX11 side:

A shader resource view that allows us to read the source texture.
A shader resource view for our constant buffer.
A large buffer to store all of the per-tile histograms. Its size will be numTiles * numHistogramBins.
An unordered access view for the per-tile histogram buffer.

We can then dispatch our compute shader. Note, that we have to round up the image size to be a multiple of our compute size; we must then be careful to only sample the real pixels. Here's the code:

#define NUM_HISTOGRAM_BINS (64)
#define HISTOGRAM_TILE_SIZE_X (32)
#define HISTOGRAM_TILE_SIZE_Y (16)

Texture2DMS sourceTextureSRV;

cbuffer toneMapGammaConstants
{
    uint outputWidth;
    uint outputHeight;

    float inputLuminanceMin;
    float inputLuminanceMax;

    uint numPerTileHistograms;

    float outputLuminanceMin;
    float outputLuminanceMax;
};

// RGB to luminance
float rgbToLuminance(float3 colour)
{
     // Use this equation: http://en.wikipedia.org/wiki/Luminance_(relative)
     return dot(float3(0.2126f, 0.7152f, 0.0722f), colour);
}

// RGB to log-luminance
float rgbToLogLuminance(float3 colour)
{
     return log(rgbToLuminance(colour));
}
// Convert a luminance to a histogram bin index
uint rgbToHistogramBin(float3 colour)
{
    // Convert to luminance
    float luminance = rgbToLogLuminance(colour);
    // Work out response curve response
    return uint(float(NUM_HISTOGRAM_BINS) * saturate((luminance - inputLuminanceMin) / (inputLuminanceMax - inputLuminanceMin)));
}

// Thread groups correspond to a tile.
// Threads correspond to a pixel

groupshared uint histogram[NUM_HISTOGRAM_BINS];

Buffer perTileHistogramSRV;
Buffer mergedHistogramSRV;
RWBuffer perTileHistogramUAV;
RWBuffer mergedHistogramUAV;

[numthreads(HISTOGRAM_TILE_SIZE_X, HISTOGRAM_TILE_SIZE_Y, 1)]
void compute_per_tile_histogram_cs(uint3 globalIdx : SV_DispatchThreadID, uint3 localIdx : SV_GroupThreadID, uint3 groupIdx : SV_GroupID)
{
    uint localIdxFlattened = localIdx.x + localIdx.y * HISTOGRAM_TILE_SIZE_X; // Which pixel within a tile
    uint tileIdxFlattened = groupIdx.x + groupIdx.y * GetNumTilesX(); // Which tile from the screen

    // Initialise the contents of the group shared memory. Only one thread does it to avoid bank conflict overhead
    if(localIdxFlattened == 0)
    {
        for(uint index = 0; index < NUM_HISTOGRAM_BINS; ++index)
        {
            histogram[index] = 0;
        }
    }

    GroupMemoryBarrierWithGroupSync();

    // For each thread, update the histogram

    // It is possible that we emit more threads than we have pixels. This is caused due to rounding up an image to a multiple of the tile size
    if(globalIdx.x < outputWidth && globalIdx.y < outputHeight)
    {
        // We just use sample 0 for efficiency reasons. It is unlikely that omitting the extra samples will make a significant difference to the histogram
        float3 pixelValue = sourceTextureSRV.Load(int2(globalIdx.xy), 0).rgb;
        uint bin = rgbToHistogramBin(pixelValue);
        InterlockedAdd(histogram[bin], 1u);
    }

    // Thread 0 outputs this thread group's histogram to the buffer
    GroupMemoryBarrierWithGroupSync();
    if(localIdxFlattened == 0)
    {
        // This could write uint4s to the UAV as an optimisation
        uint outputHistogramIndex = NUM_HISTOGRAM_BINS * tileIdxFlattened;
        for(uint index = 0; index < NUM_HISTOGRAM_BINS; ++index)
        {
            uint outputIndex = index + outputHistogramIndex;
            perTileHistogramUAV[outputIndex] = histogram[index];
        }
    }
}

We now need to take this buffer and combine it to form a single histogram for the whole image.

Merging the Per-Tile Histograms

This is a relatively simple compute shader. In this case, each thread maps to one of the tiles just computed, and each threadgroup therefore represents a group of tiles. This shader atomically adds each tile's histogram to the global histogram. In addition to the previous resources, we also now need:

A buffer to hold the histogram for the whole image.
An unordered access view to read said buffer.
A shader resource view for the per-tile histogram buffer.

// A thread per tile

#define NUM_TILES_PER_THREAD_GROUP 768 // DX10 maximum
[numthreads(NUM_TILES_PER_THREAD_GROUP, 1, 1)]
void merge_histograms_cs(uint3 globalIdx : SV_DispatchThreadID, uint3 localIdx : SV_GroupThreadID, uint3 groupIdx : SV_GroupID)
{
    // Each thread has the job of adding in the contents of the per-tile histogram to the overall histogram stored in GSM
    // This is an atomic operation, because many, many threads are trying to hit that memory simultaneously
    uint tileIndex = globalIdx.x;
    if(tileIndex < numPerTileHistograms)
    {
        for(uint binIndex = 0; binIndex < NUM_HISTOGRAM_BINS; ++binIndex)
        {
            InterlockedAdd(mergedHistogramUAV[binIndex], perTileHistogramSRV[tileIndex * NUM_HISTOGRAM_BINS + binIndex]);
        }
    }
}

Note that here we perform atomic operations on the merged histogram. You may do this on both group shared memory and unordered access views; clearly there is extra overhead in the latter. We now have a histogram for the whole image. Our job is now to convert that to a response curve. This is a function that maps from the input HDR image to a new LDR image. Wikipedia offers an easy explanation.

Computing the Response Curve

Again, more resources are needed:

A buffer holding a float for each histogram bin.
An unordered access view to write to it with.
A shader resource view for the whole-image histogram.

There is a small wrinkle in this code. We also include an extra step - the adjustment phase discussed in Page 14 of the Ward '97 paper. This is included only as an example to get you started. I use a 1x1x1 compute shader here as this is such a minimal operation that parallelism would gain little.

RWBuffer responseCurveUAV;

groupshared float frequencyPerBin[NUM_HISTOGRAM_BINS];
groupshared float initialFrequencyPerBin[NUM_HISTOGRAM_BINS];
[numthreads(1, 1, 1)]
void compute_histogram_response_curve_cs(uint3 globalIdx : SV_DispatchThreadID, uint3 localIdx : SV_GroupThreadID, uint3 groupIdx : SV_GroupID)
{
// Compute the initial frequency per-bin, and save it
float T = 0.0f;
for(uint bin = 0; bin < NUM_HISTOGRAM_BINS; ++bin)
{
float frequency = float(mergedHistogramSRV[bin]);
frequencyPerBin[bin] = frequency;
initialFrequencyPerBin[bin] = frequency;
T += frequency;
}

// Naive histogram adjustment. There are many, many such histogram modification algorithms you may seek to employ - this is
// an example implementation that will no doubt later be changed
// This is an implementation of page 14 of "A Visibility Matching Tone Reproduction Operator for High Dynamic Range Scenes"
// There are other, better approaches, like Duan 2010: http://ima.ac.uk/papers/duan2010.pdf
float recipDisplayRange = 1.0f / (log(outputLuminanceMax) - log(outputLuminanceMin));

// Histogram bin step size - in log(cd/m2)
float deltaB = (inputLuminanceMax - inputLuminanceMin) / float(NUM_HISTOGRAM_BINS); // Luminance values are already log()d

float tolerance = T * 0.025f;
float trimmings = 0.0f;
uint loops = 0;
do
{
// Work out the new histogram total
T = 0.0f;
for(uint bin2 = 0; bin2 < NUM_HISTOGRAM_BINS; ++bin2)
{
T += frequencyPerBin[bin2];
}

if(T < tolerance)
{
// This convergence is wrong - put it back to the original
T = 0.0f;
for(uint index = 0; index < NUM_HISTOGRAM_BINS; ++index)
{
frequencyPerBin[index] = initialFrequencyPerBin[index];
T += frequencyPerBin[index];
}
break;
}

// Compute the ceiling
trimmings = 0.0f;
float ceiling = T * deltaB * recipDisplayRange;
for(uint bin3 = 0; bin3 < NUM_HISTOGRAM_BINS; ++bin3)
{
if(frequencyPerBin[bin3] > ceiling)
{
trimmings += frequencyPerBin[bin3] - ceiling;
frequencyPerBin[bin3] = ceiling;
}
}
T -= trimmings;

++loops;
}
while(trimmings > tolerance && loops < 10);

// Compute the cumulative distribution function, per bin
float recipT = 1.0f / T;
float sum = 0.0f;
for(uint bin5 = 0; bin5 < NUM_HISTOGRAM_BINS; ++bin5) // By now you are wondering about bin1, bin2, bin3 etc. It is a shader compiler thing.
{
float probability = frequencyPerBin[bin5] * recipT;
sum += probability;
responseCurveUAV[bin5] = sum;
}
}

And we now have a response curve to apply to our image!

Applying the Response Curve

Again, we need one more new resource: A shader resource view to access the response curve buffer. The code is again quite simple. Here, we're running a compute shader over tiles of the image. I've also included an MSAA resolve and gamma correction step, as these are features you will probably need to implement too in a next-generation linearly correct renderer. Note, that you ideally tone-map before the MSAA resolve. I felt it was simpler to present the code this way around. Secondly, ideally your frame buffer would have separated chrominance from luminance and you would just scale the luminance. Again, I've presented it this way for a simple example.

#define TONE_MAP_GAMMA_TILE_SIZE (8)

Buffer responseCurveSRV;
RWTexture2D outputUAV;

[numthreads(TONE_MAP_GAMMA_TILE_SIZE, TONE_MAP_GAMMA_TILE_SIZE, 1)]
void main_cs(uint3 threadID : SV_DispatchThreadID)
{
float3 sourceRGB = 0;
    for(uint sampleIndex = 0; sampleIndex < NUM_MSAA_SAMPLES; ++sampleIndex)
    {
        sourceRGB += sourceTextureSRV.Load(int2(threadID.xy), sampleIndex).rgb;
    }
    sourceRGB /= NUM_MSAA_SAMPLES;

    // Work out response curve response
    uint bin = rgbToHistogramBin(sourceRGB);
    float scalar = bin < NUM_HISTOGRAM_BINS ? responseCurveSRV[bin] : 1.0f;
    sourceRGB *= scalar;

outputUAV[threadID.xy] = float4(saturate(pow(sourceRGB, 1.0f / 2.2f)), 1.0f));}

That's It

And there we have it. A simple, bare-bones implementation to illustrate how these things may be implemented. Something to get the practitioner up and running.

Distant Shadowing "brainfart"

2013-04-25T14:51:00.004+01:00

Recently, I've been having a little bit of a think about how to handle realtime distant shadowing in games, in order to avoid shadow baking.

Clearly, it's entirely possible to use some kind of cascaded shadow map extending all the way out. But, this soon ends up getting expensive in terms of textures and rendering. Alternatively, you might consider using something like stencil shadow volumes, and using the extended compute capabilities of modern GPUs to clip down the shadow volume so it doesn't consume enormous amounts of rasterisation or clipping time.

Alternatively, I've been thinking a lot about reprojection. Let's say you interrogate your cascaded shadow map and use that to build a screen-space mask of the shadow response. You want soft(ish) shadows, so you commit a lot of PCF samples.

How about reprojecting last frame's shadow mask into this frame's shadow mask, then filling in the holes with some simple filtering operations? Distant shadows are typically pretty low-frequency things, so you could well get away with some blurring and imperfection in the shadow mask. Then, you need only render the geometry that became newly visible this frame into the shadow map. You could just subtract one frame's light frustum from the other, yielding a convex polytope, and only render the intersecting objects.

Just a brainfart...

Thoughts on VPL Generation Using Cube Maps

2013-03-28T15:47:00.000+00:00

Recently, I've been doing a lot of research and thinking about realtime GI and its related ideas. A common theme in realtime GI is the use of VPLs, a way to approximate the radiance bouncing off a surface after interacting with its material. A VPL is only formed where there is a definite light-to-surface interaction. The VPL can therefore be used to simulate light bouncing off geometry.

One technique in use today for generating VPLs is the Reflective Shadow Map algorithm. Here, we render the scene from the light's point of view, and every pixel in the output image represents a surface receiving primary illumination from the light source. This is obviously a fairly efficient way of finding all the surfaces that receive radiance from a given light source without too much wastage. You can then construct VPLs using that image, reduce the working set and inject them into some lighting algorithm to bounce light around your scene.

The difficulty is, this may not scale well as each new light source would require a new RSM render.

One alternative I considered was the possibility of using a cube map to sample the primary illumination from multiple lights simultaneously. You would put your cube map at your camera position, or some other meaningful point-of-interest. You then render the scene into each cube map face. As you render it, you light it from all of your light sources simultaneously, with shadow maps too if desired. Obviously this would be an expensive render.. but it would only be done once, you can vary the resolution of it and you can always use techniques like tiled lighting to accelerate it if you have many lights.

When complete, this cube map holds the radiance at the first light bounce for many surfaces in the scene near to the point of interest - ie the camera. It contains samples for surfaces all around - both in front of and behind the camera, which is important for GI. You may then use this cube map to generate VPLs for further bounce lighting.

Now, this is as-yet a completely untested idea. I will get round to it, though. I thought I'd throw it out there as an idea for using a commonly-accelerated rendering structure to try to help the scalability of realtime GI solutions.

C++ 11 Experiences

2013-02-19T12:19:00.002+00:00

I've recently been developing a lot of code in C++ 11. I find it a huge and useful improvement over previous version of C++; here's some observations on some of the features in it that I've found rather valuable.

emplace

This simple, easily-overlooked addition to the STL has to be one of the most practically useful features of C++ 11. In itself, it does very little: it constructs a new object in-place inside a container. This in itself is very useful for optimisation purposes. However, the real advantage lies in semantics, rather than optimisation.

emplace() means that you are no longer forced to include copy semantics into classes that really, really do not want them in order to place (copy) them into a container type. If the container implements some emplace() functionality, you are free to then delete all copy assignment and copy construction operations, only providing move assignment and construction instead.

This is a huge boon for classes such as resource classes, where assignment or copy can sometimes be a very unnatural and dangerous operation. With emplace() and move semantics, you are then free to create resource classes that make extensive use of const and RAII, ensuring that only one class ever binds to a particular resource.

Lambdas for Complex Initialisation

How often have you written code such as:

SomeType blah = unwantedDefaultValue;
if(x)
{
    blah = a; // In reality this would be much more complex
}
else if(y)
{
    blah = b;
}
else
{
    blah = c;
}

Which means that "blah" is not const, when it really needs to be, and you have a big chunk of complex code polluting the flow of some code. Typically, you even have the same thing near-duplicated.

You can always factor the big initialisation statement into a function, and assign it to a const variable... but often to do so artificially externalises details you may prefer not to.

With lambdas, you can now factor all that into a lambda expression, capturing nearby variables or passing different values as a parameter:

const auto initialiser = []()
{
    if(x)
    {
        return a;
    }
    else if(y)
    {
        return b;
    }
    else
    {
        return c;
    }
};

...

const SomeType blah = initialiser();
const SomeType blah2 = initialiser2(x, y);
const SomeType blah3 = initialiser3({1, 2, 3, 4});
const SomeType blah4 = initialiser3({});
And so on. This can make the code a little more focused, more const-safe, less duplicated and easier to step through and debug.

Avoid Ternary Expressions Using Lambdas

A simple one. Ternary expressions can be a little difficult to debug, especially when they are nested:

const int a = d < e ? f :
g < h ?: i :
j < k ? m : n;

And so on. Whilst I don't advocate it, it does occur. Along the lines of the above point, using a lambda creates code you can place breakpoints on and debug easily, at no performance loss.

Reduction of Code Duplication Using Lambdas

Lambdas are very useful to reduce internal code duplication within a function. This perhaps seems a little excessive... but it can help because it factors out the details of a particular operation into one place, leaving the rest of the code simpler. Again, this avoids the problem of artifical external functions that reduce encapsulation... though I do not propose that everything be written as giant functions containing lambdas!

What's more, I've found that the compiler is generally smart enough to inline the expansion of these lambas as appropriate. So there is no performance penalty!

Move Semantics For Resource Assignment

A common issue I face with code is in creation code, where one function creates a resource, passes it to another object which is assumed to take responsibility for its lifetime management and destruction. However, often the constructor just takes a pointer to an object - this declaration is less than clear on what it will do with that resource. This doesn't tell you whether it'll just hold that reference, or manage it.

I find move semantics can help here, a little. A constructor which has an rvalue reference to some resource declares that it will acquire the resource, hold a reference to it, but also divest the calling code of that object. This perhaps makes it a little clearer that ownership of the resource is being transferred, not simply passed.

Adding the secret sauce of std::unique_ptr<> makes lifetime management even more painless.

Auto

Auto seems one of the potentially more divisive additions to C++ 11. At first, I was unsure. I can see how it can avoid template craziness... but would the lack of explicit type names make things harder to debug?

So far, I generally use auto by default as much as I possibly can. I find that if your code is making good use of functions and constructors, they tend to make it fairly clear what types are being constructed. Furthermore, it makes the code much more amenable to change and refactoring, as it reduces the number of times you have to name a particular type. This seems to be a good thing.

The main downsides of it are that you have to be extremely careful to remember to include & to make some types a reference. If you do not, then you end up implicitly copying objects around, which can lead to performance, correctness or compilation issues.

Range-based for

This one's not such a big deal. It can be useful, as it can save on code verbosity... but so often, you need to iterate over multiple containers simultaneously, and that's where range-based for ceases to be useful.

constexpr

This is a nice simple win. This at least means you don't have to start making trivial functions into templates to get some compile-time optimisations. It's unfortunate that it insists on the single-line return statements, as this can mean things such as data tables spill out of the scope of the function.

Class Member Initialisation

I find that being able to assign values to class instance members inside the class declaration to be enormously useful. This allows you to vastly simplify your constructors, especially when you have multiple constructors. This removes a lot of duplication, and allows you to instead implement only the specialisations. This seems to be the right way round to me.

Promises, and Futures

By and large, these are quite nice facilities. The big thing they lack though is the ability to poll a future. This seems like quite a major omission, and can make them difficult to use for some algorithms. As such, I tend to find myself rolling my own cut-down variants for many purposes.

In summary, it's very agreeable. No real change is mandated, so you simply opt in or opt out of as many or as few new features you'd like to employ. I'm a convert.

Job Graphs in C++ 11

2012-11-02T21:30:00.005+00:00

Introduction

I've recently been experimenting with a lot of the new features in C++ 11, purely for discovery and experimentation. One major new feature are the threading libraries. As using is the best way of learning, I set upon the task of writing a currently-fashionable task/job-graph system using these new features.

I'm not promoting this algorithm as the one true approach for parallelism, as there are many other valuable algorithms that are much more suitable in certain occasions. This one seems a useful, generic solution that fits many cases fairly well. In the process, I found a few interesting quirks during implementation that I thought may be interesting for others.

Job Graphs

Job graphs are currently finding wide favour as an approach to the challenge of building a flexible architecture for multithreaded processing. In this scheme, jobs are stored in a graph structure. The edges of the graph are used to represent ordering and data dependencies between jobs. Assuming the graph is structured correctly, this ensures that work is processed in the correct order to ensure the correct and timely flow of data. Nodes that do not share a data dependency may be executed concurrently. The programmer shoulders the responsibility of structuring the graph correctly, and ensuring any conflicts between resources are either resolved via the graph, or using some locking mechanism.

Typically, the jobs are pushed into some kind of queue and distributed to a thread pool. Ideally you would also use work-stealing for maximum efficiency.

The goal is to use this graph structure to allow non-dependent work to execute concurrently, saturating all hardware threads. It is a very simplistic idea, yet at its most it is NP-hard to solve.

C++ 11 Threading Overview

The C++ 11 threading library is really quite a lightweight, spare API. Broadly, it offers the facilities to:

Create and manage threads in a cross platform fashion.
Synchronise activity between these threads using simple primitives such as mutexes, condition variables and promises/futures.
Package up tasks to be executed on a thread.

Implementation

The implementation of a job graph requires relatively few components:

A job class
A job graph container
A job pool
A thread pool

Job Class

The job class is fairly straightforward. It serves two purposes: to package up some functionality to be executed, and have knowledge of precursor jobs and dependent jobs (this concept is sometimes described as continuations).

The major consideration in the design of the job class is the mechanism used to sequence their operation. Ideally, you want to use some kind of OS threading primitive as opposed to some kind of user-space mechanism. The OS is in a far better position to manage threads as it has clear information of the state and dependencies of all of them, and can wake and sleep threads with minimal context switches. User space code often ends up with some kind of inefficient wake-predicate-sleep polling pattern causing a lot of wasted time and context switches.

I therefore chose to use a C++ 11 thread library primitive: std::promise and std::future. A promise provides the ability to calculate some value asynchronously, in a deferred context, in the future. It represents a contract for some piece of code to provide a value based on some inputs. A future is the mechanism used to retrieve that value. The promise is the producer, the future is the client.

The promise and future abstraction provides a good mechanism to couple jobs together. A job can provide a "future" to each of its dependent jobs. The dependent job uses this future to wait for the job's completion. The future could pass data, or it could solely be used to sequence computation.

A second feature of C++ 11 that proves useful for the job class is std::function and lambda functions. An std::function is simply a wrapper for some callable target. My job class contains one of these to hold the job's processing kernel. It is initialised in the constructor, which allows me to pass in a lambda function when creating the job. This provides an interesting alternative to using inheritance and virtual functions, and the resultant proliferation of classes. Now, I have just one Job class, and the data handles the variation.

Job.h:

class Job;

//! This monitors the completion of a number of jobs, and sends out a signal when they're all done

class JobCompletionMonitor

{

public:

JobCompletionMonitor(void) :

m_numJobsPending(0)

{

}

// Forbidden constructors

JobCompletionMonitor(JobCompletionMonitor&) = delete;

JobCompletionMonitor& operator=(JobCompletionMonitor&) = delete;

// Jobs tell the monitor what's going on

void notifyStart(Job* job);

void notifyFinish(Job* job);

// Other people can find out when they're all done

void waitForCompletion(void);

bool allJobsFinished(void) const;

private:

std::atomic_int m_numJobsPending;

std::condition_variable m_allJobsComplete;

std::mutex m_mutex;

};

//! A self-contained job that is executed somewhere in the frame

class Job : public objects::CoreObject

{

public:

// Priority system

typedef unsigned int Priority;

static const Priority HighPriority = 1000;

static const Priority AboveNormalPriority = 5000;

static const Priority NormalPriority = 10000;

static const Priority BelowNormalPriority = 15000;

static const Priority LowPriority = 20000;

// Permitted constructors

Job(const resources::NameManager::Name& name,

const std::initializer_list<std::shared_ptr<Job>>& precursorJobs = {},

std::function<void(float)> kernel = [](const float){},

const Priority priority = NormalPriority);

~Job();

// Forboidden constructors

Job(Job& job) = delete;

Job& operator=(Job& job) = delete;

// Management of dependents

void addDependent(const std::shared_ptr<Job>& dependent);

void removeDependent(const std::shared_ptr<Job>& dependent);

// Execution API

void prepare(JobCompletionMonitor* monitor);

bool canRun(void) const;

void execute(const float timeStep);

void cleanup(void);

// Accessors

inline const std::vector<std::shared_ptr<Job>>& getDependents(void) const

{

return m_dependents;

}

inline std::vector<std::shared_ptr<Job>> getPrecursorJobs(void) const

{

return m_precursors;

}

inline bool operator<(const Job& other) const

{

return m_priority < other.m_priority;

}

private:

// Connection to other jobs

std::vector<std::shared_ptr<Job>> m_precursors; //!< The job that must complete before we can run. We are a continuation of this job

std::vector<std::shared_ptr<Job>> m_dependents; //!< A list of all the jobs that continue after us

// Futures to sequence processing

std::vector<std::future<int>> m_precursorFutures; //!< The future that is used to tell us when to start - if we have a precursor

std::vector<std::promise<int>> m_dependentPromises; //!< The promises we fulfill to trigger other continuations

// Override this function to implement your logic

std::function<void(const float)> m_kernel; //!< A simple routine to call to implement the logic for this task

// What is the execution priority of this job?

Priority m_priority; //!< Execution priority of this job

JobCompletionMonitor* m_completionMonitor = nullptr; //!< This object monitors the job. The job must notify it on completion

};

Job.cpp:

//! Have all of our jobs now finished?

bool JobCompletionMonitor::allJobsFinished(void) const

{

return m_numJobsPending == 0;

}

//! The job tells the monitor when it starts

void JobCompletionMonitor::notifyStart(Job* job)

{

INC_RT_ASSERT(job != nullptr);

++m_numJobsPending;

}

//! The job notifies the monitor that it is complete

void JobCompletionMonitor::notifyFinish(Job* job)

{

INC_RT_ASSERT(job != nullptr);

--m_numJobsPending;

if(m_numJobsPending == 0)

{

m_allJobsComplete.notify_all();

}

//! Wait until all the jobs managed by this monitor are complete

void JobCompletionMonitor::waitForCompletion(void)

{

std::unique_lock<std::mutex> lock(m_mutex);

m_allJobsComplete.wait(lock);

INC_RT_ASSERT(m_numJobsPending == 0);

}

//! Construction

Job::Job(const resources::NameManager::Name& name,

const std::initializer_list<std::shared_ptr<Job>>& precursorJobs,

std::function<void(float)> kernel,

const Priority priority) :

objects::CoreObject(name),

m_precursors(precursorJobs),

m_dependents(),

m_dependentPromises(),

m_kernel(kernel),

m_priority(priority)

{

}

//! Destruction

Job::~Job()

{

// Add some assertions here to ensure the job has been correctly deleted

}

// Notify this job of a new dependent that expects to be told when this has finished

void Job::addDependent(const std::shared_ptr<Job>& dependent)

{

INC_RT_ASSERT(dependent != nullptr);

// If this job isn't already a dependent, then add it

if(std::find(m_dependents.begin(), m_dependents.end(), dependent) == m_dependents.end())

{

m_dependents.push_back(dependent);

}

// Remove a dependency from this job

void Job::removeDependent(const std::shared_ptr<Job>& dependent)

{

INC_RT_ASSERT(dependent != nullptr);

// If this job isn't already a dependent, then add it, and re-issue futures

auto it = std::find(m_dependents.begin(), m_dependents.end(), dependent);

if(it != m_dependents.end())

{

// Remove this dependent

m_dependents.erase(it);

}

//! Reset the job ready for graph traversal and execution

void Job::prepare(JobCompletionMonitor* monitor)

{

INC_RT_ASSERT(monitor != nullptr);

m_completionMonitor = monitor;

// Ask our precursors for a future

for(auto& precursor : m_precursors)

{

precursor->m_dependentPromises.push_back(std::promise<int>());

m_precursorFutures.push_back(precursor->m_dependentPromises.back().get_future());

}

//! Clean up after execution. Remove all the promises and futures

void Job::cleanup(void)

{

m_precursorFutures.clear();

m_dependentPromises.clear();

}

// Can this job run? Have all of our futures been satisfied?

bool Job::canRun(void) const

{

for(auto& it : m_precursorFutures)

{

if(it.valid() == false)

{

return false;

}

return true;

}

// Actually run this job

void Job::execute(const float timeStep)

{

// Wait to receive the promised values from all precursor jobs

for(auto& future : m_precursorFutures)

{

future.wait();

(void)future.get();

}

// Do our actual work

m_kernel(timeStep);

// Fulfill promises to our children

for(auto& promise : m_dependentPromises)

{

promise.set_value(1);

}

m_completionMonitor->notifyFinish(this);

// Now that we are done, clear out all of the stale old futures we were waiting for -

// prepare() will re-issue them next time around

m_precursorFutures.clear();

}

Job Graph

The job graph really has few responsibilities: To contain the jobs, traverse the graph and manage the resources needed for execution.

The job graphs are contained in a simple graph structure. I may later shift to a more capable container such as boost::graph.

The main task at hand for the job graph is to traverse the graph, pushing all of the jobs into an ordered pool ready for execution. A thread pool then pulls jobs from this pool when they are ready for execution.

JobGraph.h:

// A graph structure holding all the jobs and conditions

class JobGraph

{

public:

JobGraph(const size_t numThreads);

~JobGraph();

void add(const std::shared_ptr<Job>& job);

void remove(const std::shared_ptr<Job>& job);

void execute(const float timeStep);

private:

// At the start of the frame, the job graph is traversed and all the jobs to be processed are placed into this pool

JobPool m_jobPool;

// This is a pool of threads which work on executing those jobs

// Threads are free to steal jobs from the overall pool to keep things ticking forwards

ThreadPool m_threadPool;

// This is a special root-level job

std::shared_ptr<Job> m_rootJob;

// Need to measure when all jobs have completed

// (not just dispatched - *completed*)

JobCompletionMonitor m_jobCompletionMonitor;

// Mutually exclude operations like construction, adding jobs and execution

std::mutex m_mutex;

// Does a (parent) job exist within the graph?

std::unordered_map<Job*, bool> m_jobContained;

};

JobGraph.cpp:

// A generic minimal visitor pattern

class JobGraphVisitor

{

public:

virtual ~JobGraphVisitor() {}

void traverse(Job* job)

{

INC_RT_ASSERT(job != nullptr);

visit(job);

for(auto child : job->getDependents())

{

traverse(child.get());

}

};

private:

virtual void visit(Job* job) {}

};

// Visit with a view to pushing a job to the task pool

class PushJobVisitor : public JobGraphVisitor

{

public:

inline PushJobVisitor(JobPool& jobPool, const float timeStep, JobCompletionMonitor* jobCompletionMonitor) :

m_jobPool(jobPool),

m_timeStep(timeStep),

m_jobCompletionMonitor(jobCompletionMonitor)

{

}

private:

JobPool& m_jobPool;

const float m_timeStep;

JobCompletionMonitor* const m_jobCompletionMonitor;

std::unordered_map<Job*, bool> m_jobPushed;

virtual void visit(Job* job) override

{

if(m_jobPushed[job] == false)

{

m_jobCompletionMonitor->notifyStart(job);

m_jobPool.push(job);

m_jobPushed[job] = true;

}

};

// Reset jobs in preparation

class ResetJobVisitor : public JobGraphVisitor

{

public:

inline ResetJobVisitor(JobCompletionMonitor* jobCompletionMonitor) :

m_jobCompletionMonitor(jobCompletionMonitor)

{

}

private:

JobCompletionMonitor* const m_jobCompletionMonitor;

virtual void visit(Job* job) override

{

job->prepare(m_jobCompletionMonitor);

}

};

//! Post-execution by-product clean-up

class CleanUpJobVisitor : public JobGraphVisitor

{

public:

inline CleanUpJobVisitor(void)

{

}

private:

virtual void visit(Job* job) override

{

job->cleanup();

}

};

// Construction

JobGraph::JobGraph(const size_t numThreads) :

m_jobPool(),

m_threadPool(numThreads, m_jobPool),

m_rootJob(new Job(resources::getUniqueName("Root Job")))

{

INC_RT_ASSERT(numThreads != 0);

std::unique_lock<std::mutex> lock(m_mutex);

}

// Destruction

JobGraph::~JobGraph()

{

INC_RT_ASSERT(m_jobPool.jobAvailable() == false);

m_jobPool.shutdown();

m_threadPool.shutdown();

std::this_thread::yield();

}

// Add a job to the graph

void JobGraph::add(const std::shared_ptr<Job>& job)

{

INC_RT_ASSERT(job != nullptr);

std::unique_lock<std::mutex> lock(m_mutex);

if(job->getPrecursorJobs().size() == 0)

{

// A null precursor means that it's a root level task that lives under our

// special root job

m_rootJob->addDependent(job);

}

else

{

// Assert to ensure that the parent job actually exists within the graph!

for(auto& precursor : job->getPrecursorJobs())

{

INC_RT_ASSERT(m_jobContained[precursor.get()]);

precursor->addDependent(job);

}

m_jobContained[job.get()] = true;

}

// Remove a job from the graph

void JobGraph::remove(const std::shared_ptr<Job>& job)

{

INC_RT_ASSERT(job != nullptr);

std::unique_lock<std::mutex> lock(m_mutex);

if(job->getPrecursorJobs().size() == 0)

{

m_rootJob->removeDependent(job);

}

else

{

for(auto& precursor : job->getPrecursorJobs())

{

precursor->removeDependent(job);

}

m_jobContained[job.get()] = false;

}

// Execute all jobs in the graph

void JobGraph::execute(const float timeStep)

{

std::unique_lock<std::mutex> lock(m_mutex);

// Make sure the job queue is fully empty from last time

INC_RT_ASSERT(m_jobPool.jobAvailable() == false);

// Get it all reset ready for execution

{

ResetJobVisitor resetJobVisitor(&m_jobCompletionMonitor);

resetJobVisitor.traverse(m_rootJob.get());

}

// Traverse the whole job graph, pushing the jobs into the task pool

{

PushJobVisitor pushJobVisitor(m_jobPool, timeStep, &m_jobCompletionMonitor);

pushJobVisitor.traverse(m_rootJob.get());

}

// (Note that tasks will be consumed asynchronously as soon as we start to push them)

m_jobCompletionMonitor.waitForCompletion();

// We should be all done now. Nothing should be here

INC_RT_ASSERT(m_jobPool.jobAvailable() == false);

// Get it all cleaned up ready for next time

{

CleanUpJobVisitor cleanUpVisitor;

cleanUpVisitor.traverse(m_rootJob.get());

}

Job Pool

It is initially tempting to simply traverse the graph, push all of the jobs out for execution, and let the OS arbitrate everything. There are a few difficulties here. First of all, unless your target platform's runtime explicitly models tasks or jobs, this tends to lead to mass thread creation and destruction - a substantial overhead. Secondly, you may wish to cancel the execution of jobs after a precursor has completed. Thirdly, it makes little sense to dispatch many jobs for execution that have no chance of execution. They will just immediately hit a threading hazard and context-switch away. Finally, consider the case where you have N cores and N+1 jobs. Job 0 does some work, and jobs 1..N wait on job 0. Imagine if your target platform has a pooly-implemented runtime, and assigns jobs 1..N to cores instead of job 0. Those jobs may wait, and job 0 may not execute. This case is unlikely in a well-implemented runtime, but it is a risk, and why present code that is inherently at risk?

The job pool is simply an std::set of jobs, ordered by job priority. Whenever a thread in the thread pool wakes up, it pops the first (highest priority) job from the pool and executes it.

JobPool.h:

// This represents a pool of tasks that need doing, for all threads to

// pull from (and therefore work-steal from)

class JobPool

{

public:

JobPool(void);

JobPool(JobPool&) = delete;

JobPool& operator=(JobPool&) = delete;

bool jobAvailable(void);

void push(Job* job);

Job* pop(void);

void waitForJob(void);

bool shouldTerminate(void);

void shutdown(void);

private:

//! This gets notified when a task is pushed, to allow the other consuming

//! threads to wake up

std::condition_variable m_jobPushNotification;

std::condition_variable m_queueEmptyNotification;

//! We maintain a set of jobs, so that we can find the first, highest-priority job

//! that is in an executable state

std::set<Job*> m_jobs;

//! This mutex is used to ensure thread safety when pushing a job in or pulling it out

std::mutex m_mutex;

//! Used when it's time to shut down

std::atomic<bool> m_shouldTerminate;

};

JobPool.cpp:

// Construct it

JobPool::JobPool() :

m_shouldTerminate(false)

{

}

// A consumer calls this function to wait for notification via the condition_variable of a new task

void JobPool::waitForJob(void)

{

std::unique_lock<std::mutex> lock(m_mutex);

m_jobPushNotification.wait_for(lock, std::chrono::milliseconds(100));

// The above timeout means that spurious wakeup is possible

}

// Is there a job available in the set?

bool JobPool::jobAvailable(void)

{

std::lock_guard<std::mutex> lock(m_mutex);

return m_jobs.empty() == false;

}

// Push a new job into the pool

void JobPool::push(Job* job)

{

INC_RT_ASSERT(job != nullptr);

std::lock_guard<std::mutex> lock(m_mutex);

INC_RT_ASSERT(m_shouldTerminate == false);

m_jobs.insert(job);

m_jobPushNotification.notify_one();

}

// Try to return an executable job from the pool

Job* JobPool::pop(void)

{

std::lock_guard<std::mutex> lock(m_mutex);

if(m_jobs.empty())

{

// Nothing to run!

return nullptr;

}

else

{

// Find the highest priority job that can run

for(auto job : m_jobs)

{

if(job->canRun())

{

m_jobs.erase(job);

return job;

}

// It is possible that we're asked for a job but none are able to run... yet

return nullptr;

}

//! Should we terminate now?

bool JobPool::shouldTerminate(void)

{

return m_shouldTerminate;

}

//! Send termination message

void JobPool::shutdown(void)

{

std::lock_guard<std::mutex> lock(m_mutex);

m_shouldTerminate = true;

m_jobs.clear();

m_jobPushNotification.notify_all();

}

Thread Pool

The purpose of the thread pool is to provide a persistent set of threads capable of executing jobs, saving on the cost of thread creation and destruction. These threads pull jobs from the job pool until they are asked to terminate.

Spurious Wakeup

The threads in the thread pool wait on an std::condition_variable to notify them on the availability of a job. It is important to consider the possibility of spurious wakeup in this system. Therefore, when a thread wakes up, it must ensure a job actually exists in the pool to be executed.

I actually deliberately allow the std::condition_variable::wait() to time out after 100 ms. The rationale for this is that some operating systems become suspicious of threads that sleep for long periods of time. For stability reasons, I prefer my system to wake up unneccesarily and to service the outer thread loop than wait for infinite amounts of time. This therefore deliberately introduces some degree of spurious wakeup into the system.

ThreadPool.h:

// A pool of threads hungry to execute tasks

class ThreadPool

{

public:

ThreadPool(const size_t numThreads, JobPool& jobPool);

~ThreadPool();

void shutdown(void);

private:

//! The set of threads in the thread pool

std::vector<std::thread> m_threads;

};

ThreadPool.cpp:

// Construction / destruction

ThreadPool::ThreadPool(const size_t numThreads, JobPool& jobPool)

{

for(size_t thread = 0; thread < numThreads; ++thread)

{

m_threads.push_back(std::thread([](JobPool* jobPool)

{

for(;;)

{

// Sleep, until a job is available

// (use a condition_variable here)

jobPool->waitForJob();

// Have we been told to shut down?

if(jobPool->shouldTerminate())

{

return;

}

// Pull a job out of the task pool and execute it

if(jobPool->jobAvailable())

{

// Whilst we'd prefer to execute task() here, instead we have to do

// the pop-and-execute as an atomic-like operation to avoid

// two threads popping the same task, or one thread popping

// a task that the other thread expected to receive

Job* job = jobPool->pop();

if(job != nullptr)

{

// Clearly the correct time step should be plumbed here

job->execute(1.0f / 30.0f);

}

}, &jobPool));

}

ThreadPool::~ThreadPool()

{

}

// Shutdown the threads. You must have told the JobPool to terminate

void ThreadPool::shutdown(void)

{

for(auto& thread : m_threads)

{

thread.join();

}

Findings

The new C++ 11 libraries largely made this system easy and quick to implement. The libraries are quite clearly defined, unobtrusive, easy to use with few dependencies.

One of the main difficulties encountered in the implementation of this system is simply the life cycles of the C++ 11 std::promises and std::futures. Their life cycles are not clearly documented and they can be difficult to reason about. It was only by reading the STL code that I found these primitives fundamentally hold a "use and discard" policy. This leads to the need for some code to create and destroy these resources each frame. Note that this may possibly vary across platforms, but this seems the best least-common-denominator assumption. This is largely the reason I avoided the use of std::packaged_task. With things like move semantics transferring the object between threads and queues, it becomes very difficult to reason about the life cycle of the object! Furthermore, once you have your own job pool and thread pool, it largely obviates the need for this class. This class only seems useful for "go and load some assets on a thread" type work.

The system itself has thus far proven stable. I have unit tests for a variety of graph configurations, which all pass, and it successfully manages to run my current engine. This engine is in its infancy, and so has not massively load tested this system... but, ultimately, this is meant to be an illustrative, exploratory post!

The Coder-Trap of Hypotheses

2012-05-17T21:47:00.000+01:00

As my project is getting towards the end, again, I'm doing a fair bit of debugging. Every so often, a horrible bug crops up that takes a few days, a week, two weeks to fix. Every time I have one, I try to think, on reflection, what could I have done to fix it quicker?

I've noticed a pattern emerging.

Briefly, bugs take me longer to fix when I allow myself to be seduced into chasing half-formed hypotheses as to what the bug is. These hypotheses often prove ill-founded, leading to wasted time that doesn't often add real data or evidence to your bug search. Really, those intermediate hypotheses are little more than speculation or conjecture.

I'm becoming increasingly aware that the more effective strategy is not to consider hypotheses until as late as possible; instead, concentrate on gathering data, evidence and devising test fixtures. This is a very difficult discipline as people naturally love to formulate hypotheses to explain things. They are seductive. They feel like you've found a short-cut to the solution, you've caught your quarry, the tension of the bug fix is over. You get to announce your triumph.

But it's rare that you fix the most difficult bugs with a single brilliant moment of insight. More often, it takes hours of analysis and evidence gathering. That's the most reliable way to reach a well-founded hypothesis that completely fixes the bug, rather than a partial fix, a temporary fix or a heisenbug fix.

Interruption Handling Techniques for Programmers

2012-04-12T13:53:00.001+01:00

Once more, I find myself in the post-alpha run to the finish line. And once again, I find myself getting busier and busier as the completion date looms.

One of the most challenging issues at this time is managing the stream of traffic to your desk. Other people are very busy too. People have questions. People need answers. And answering those questions takes time away from handling your own workload, allowing it to multiply like bacteria. I find programmers tend to employ several archetypal techniques to handle those interruptions:

1. Just help them. Indulge them. Offer them as much help as you can. Get their problem fixed, as it helps the team move forward. Really, this is the ideal strategy for the whole project. Someone has come to you with a problem, they may be at their wits end, you should help them fix it.

The problem is that it's a win-win-lose proposition. It's a win for the individual, it's a win for the team, but it's a lose for you.

2. Say you don't know. Strangely, this is a rare tactic. People often feel it's preferable to have a stab, and speculate as to the cause or the solution, even without any useful knowledge or data. Email is rife with this kind of thing. This isn't actually useful. The individual with the problem is likely to have conducted much the same process and reached similar conclusions. It's nice that you're trying, but this just consumes time and doesn't advance a solution. Better to say you don't know, and refer them to someone else who may be able to help.

3. Pretend you don't know. You may know the answer, but since fully helping is a time sink for you, you may judge it more profitable to deny all knowledge and continue work.

Until someone rats you out.

Not a great plan. Not professional. Instead, it's better to admit knowledge where you have it, but if you have a more pressing matter, explain yourself and offer help when you're done. If it's urgent, they can find help elsewhere, if not they can wait. Chances are they'll solve it anyway in the meantime, increasing their own knowledge and learning.

4. Measured help. Give them just enough help to get them moving forward again, but stop short of full disclosure. The advantages to this technique are that you're offering help, the individual is employing their own initiative and learning and everyone gets to move forward. The downside is that they may return to your desk a further 38 times before the problem is solved.

5. Be brusque. Make people not want to ask. This is a win for you, in the short term, but a lose for the team and a lose for yourself in the long term as people will not offer help in return. It can be difficult to resist this strategy when you're stressed!

6. Refer them to a subordinate. The joys of seniority! Pass the problem down the line to the next guy. It's a win for you, but it has mixed results. It may result in a win for the delegate, as they learn something and increase their standing within the team. It may be a win for the individual. It's a win for you, as it frees your time up. The problem is if the delegate cannot solve that problem. The delegate burns a lot of time, the individual with the problem loses time, and ultimately you end up having to help out anyway. Choose carefully.

7. Dazzle them with science. This is a strange one. When faced with a non-technical user, they blaze them with incredible technical knowledge and prowess. A scintillating analysis, definition and solution to the problem is offered; the individual looks back, blankly. No solution is made. The individual retreats, and contacts your line manager to request that you implement the solution instead. You clearly know what you're talking about.

8. Debate them on it. One for the orators. Any problem can be negated with a gallant argument, excellent critical analysis and undeniable logic. Your inquirer retreats, defeated. No more problem you think! No. Instead, the inquirer learns to cut you out of the loop, and ratify and implement a solution elsewhere, elsehow.

9. Blame the build server / Visual Studio / Perforce / Whatever. Clearly the programmer doesn't honestly believe this is the case. It's just a play to buy some time, or to momentarily offload the stress.

10. Counter-problem! You come to me with a problem? You're gonna leave with a bigger one! Thats the Chicago way.

R & D - It keeps you honest!

2012-02-11T20:53:00.001+00:00

I've recently been working in an R&D group, developing various features that are quite experimental and wouldn't normally be scheduled as part of a project. This is intended to be "practical" R&D - we aim to develop features that have a chance of making it into the game, as opposed to features that are blue sky and would never actually be shipped.

This kind of development approach subtly shifts your approach to developing features. Briefly, you become very focused on doing what is necessary to get the feature into the game. You focus on making sure the feature is on-budget. You make sure that you consult with affected parties. You work out the consequences of incorporating the system and seek to mitigate any problems. You suddenly find yourself becoming a salesman, trying to establish a coalition of support for the feature. Your feature really has to prove itself worthwhile so that management will entertain the perceived risk of incorporating it.

Strangely, you don't always find yourself doing these things during the development of "standard" features. Because they are planned as part of the development effort, their inclusion in the title is automatically assured - pending no major disasters. You don't have to sell them. These features don't need to earn their place. Consequently, it can be very easy to find yourself "going down a rabbit hole" with features like this - exploring functionality that is interesting, but not necessarily focused on shipping. Perhaps these features are ultimately more risk-intensive, as they won't be quite so critically examined as the R&D features.

I tend to find that the whole "posture" you adopt during the development of an R&D feature that is not automatically included in the game is ultimately much better from the perspective of careful software development. It does at least feel more like the kind of process intended by the creators of various software development methodologies. It feels very much along the lines of the approach advocated by Michael Saladino in his "Why Red Lining Your Team Leads to Ineffectiveness".

Perhaps it's not for everyone. Perhaps it's not for every team. Perhaps some may find it quite a pressured, stressful kind of approach. But I really quite like it. It certainly feels a much more considered, deliberate approach. It seems that developing features in branches, that have to prove their worth for integration, sets in motion a slightly changed, more focused set of processes than may be typical.

Sparse Voxel Octrees in Haskell

2012-02-09T13:56:00.002+00:00

It's been a long time coming (daughters birth impedes Haskell coding!), but I've finally implemented and uploaded support for Sparse Voxel Octrees in my Haskell Renderer.

It's quite an interesting algorithm to implement in something like Haskell. It's interesting because it really questions the balance of strict vs lazy evaluation in Haskell. Which is better - strict construction and strict traversal, with the platform optimisations that may bring, or lazy construction-as-you-traverse?

In an imperative language, you'd run a complete process to construct the tree, and then traverse the completed tree for each intersection query. You have to pay the total construction cost plus the traversal cost. Since you'd be constructing and storing the whole tree, you need to invest significant effort in data compression. On the other hand, the direct data representation would permit you to make a number of optimisations aimed at evaluation efficiency.

With Haskell, the lazy evaluation thunks can do some of that work for you. Instead of trying to crunch the whole tree down to be as small as possible, you may choose instead to rely on the Haskell runtime to just-in-time construct the relevant portions of it. You would only pay the cost of construction for the parts of the tree you traverse. Would this sum be lower than a strict construction and traversal?

Alternatively, In Haskell, you can still go for the strict option and implement it much as you would with C.

What's interesting is that it introduces this challenging, general question of what offers the most efficiency: A domain-specific, efficiently implemented clever algorithm, or a somewhat more generic algorithm, less specifically optimised but instead relying on lazy evaluation and the higher-order optimisations that may automatically permit?

I'm not sure what the answer is to this question right now. I suspect there is a cross over point that is very specific to the data set in question and therefore your particular application.

Code

My current Haskell SVO code is on github, in my Haskell Renderer:

https://github.com/TomHammersley/HaskellRenderer

I'm going to fix up the various remaining TODOs and then pursue a more optimised, C-like version along the lines of the various SIMD-like SVO implementations that have recently emerged.

Implementation Notes

This was ultimately quite an elegant algorithm to implement in Haskell. Much of the algorithm is built off simple geometric tests. These are quite compact and elegant to represent in a functional language like Haskell. For instance, a ray-box test is simply:

boundingBoxIntersectRay :: AABB -> Ray -> Maybe (Double, Double)
boundingBoxIntersectRay (bounds0, bounds1) (Ray rayOrg _ invRayDir rayLen)
     | tmax > tmin = Nothing
     | tmin > 0 && tmin < rayLen = Just (tmin, tmax `Prelude.min` rayLen)
     | tmax > 0 && tmax < rayLen = Just (tmax, tmax)
     | otherwise = Nothing
where
     (tmin, tmax) = foldr1 (\(a0, b0) (a1, b1) -> (a0 `Prelude.max` a1, b0 `Prelude.min` b1)) (map intercepts [vecX, vecY, vecZ])

intercepts f = let x0 = (f bounds0 - f rayOrg) * f invRayDir
                   x1 = (f bounds1 - f rayOrg) * f invRayDir
               in (x0 `Prelude.min` x1, x0 `Prelude.max` x1)

Here we simply define an operation to calculate the slab intercepts of the box, and then repeat that operation over X, Y and Z, then fold to reduce the results. We're offering a description of the algorithm to Haskell in high level primitives, permitting optimisation and reduction - as opposed to clumsily trying to indirectly tell a C++ compiler what you'd really like it to do. It's surprising how often this sort of thing proves faster in Haskell. For example, using the "any" or "all" functions has proven faster than a set of || or && clauses for me in the past. I could really imagine that a thunk-less version of Haskell, free of the irritation of considering C++ aliasing might produce some very efficient code.

The traversal code is a simple recursive traversal, with few (as yet) clever tricks. It terminates either at the leaves, at the maximum traversal depth or when the feature becomes too small when projected to screen space (though this calculation needs further elaboration). It therefore gives a simple, automatic LOD control.

And so, finally, here it is, our SVO lego-sphere:

Peak Abstraction

2012-02-02T16:46:00.003+00:00

An interesting phenomena I've observed over my years as a programmer is the occurrence of "Peak Abstraction" within a young programmer - and their subsequent improvement.

There is a common pattern. It typically occurs after 3-4 years of professional practice. A programmer is growing in his confidence. He begins to tackle more complex, challenging problems. And, he reflects those problems with complex, challenging data structures. He has not yet learnt the fine art of simplicity.

Every member of every data structure is a pointer or a reference. There are no simple data types. Everything becomes a "value dictionary", a "graph-based visitor" or a "multi-dimensional ring buffer". Disassembly of his code reveals it to be little more than an infinite sequence of lwz, lwz, lwz, lwz and more lwz. Perhaps he/she develops a hobbyist interest in Java. They may develop an editor that connects boxes with lines - such editors are great for filling with lots and lots of abstraction! Often, she/he may favour Boost; perhaps he even becomes your office's Boost Guy. Even simple plain old data types such as integers are suddenly held within a smart-pointer container and reference counted. And with a pimpl to retrieve their value.

Those around them become increasingly desperate. They can see the code becoming inflated, inefficient, unmaintainable. Spending time reading this individual's code is like spending time exposed to a radioactive isotope. You can only tolerate a few moments before exhibiting signs of sickness. You attempt to modify their code. 478 pages of template errors dissuade you from considering this foolish course of action again.

And then it happens. The programmer realises the error of their ways. They realise that they are adding neither design value nor computational value. They realise they are adding overhead and maintenance trouble. Perhaps they undergo their first encounter with profiling software. Said software may disclose the full terror of the D$ dirty bomb they have unleashed upon the code. Perhaps they are "introduced" to the profiler by an irate lead in some sort of shotgun-wedding affair. Whatever, the programmer realises that solving complex problems does not require needlessly complex data structures. Such structures impress nobody. People are impressed by simple, elegant solutions. And thus the recovery begins. Data structures become simple, expressive, and maintainable. Data cache misses dip under 20ms/frame for the first time. Co-workers enjoy modifying their code. The compiler begins to emit instructions such as "add" or "mul". Sadly the compiler may be too fond of outputting "jump" instructions at this stage... The programmer themself becomes wary of others who are pre-Peak Abstraction.

One even wonders if, like it's namesake "Peak Oil", there may be a global Peak Abstraction. Perhaps the practising body of programmers may all undergo Peak Abstraction in the next 5-10 years and we will all left never having to refactor some template monstrosity.

But then, where's the fun in that? Don't programmers just love to be indignant about others' code?

Haskell Photon Mapper Released to Hackage

2011-09-13T21:19:00.000+01:00

I've now uploaded my Haskell photon mapper and raytracer to hackage. It is listed as a package called "crocodile" (named after my son's current obsession with crocodiles!).

It's all working on the Cornell box. It's pretty memory hungry, and has a few self-declared things to fix (eg irradiance gradients), but it's out there to play with.

Enjoy.

This is what happens

2011-08-25T09:42:00.003+01:00

When universities teach Java!

(some random bit of disassembled code found on my travels)

00A5CA74 EC7E1824 fdivs f3,f30,f3                          56 PIPE
00A5CA78 FC406850 fneg f2,f13
00A5CA7C D0660010 stfs f3,0x10(r6)
00A5CA80 D0C90010 stfs f6,0x10(r9)
00A5CA84 D981FF60 stfd f12,-0xA0(r1)                        PIPE
00A5CA88 E881FF60 ld r4,-0xA0(r1)                          70 (00A5CA84) LHS[01]
00A5CA8C D961FF60 stfd f11,-0xA0(r1)                        PIPE
00A5CA90 EBA1FF60 ld r29,-0xA0(r1)                         70 (00A5CA8C) LHS[01]
00A5CA94 D941FF60 stfd f10,-0xA0(r1)                        PIPE
00A5CA98 EB81FF60 ld r28,-0xA0(r1)                         70 (00A5CA94) LHS[01]
00A5CA9C D921FF60 stfd f9,-0xA0(r1)                         PIPE
00A5CAA0 EB61FF60 ld r27,-0xA0(r1)                         70 (00A5CA9C) LHS[01]
00A5CAA4 D901FF60 stfd f8,-0xA0(r1)                         PIPE
00A5CAA8 EB41FF60 ld r26,-0xA0(r1)                         70 (00A5CAA4) LHS[01]
00A5CAAC D8E1FF60 stfd f7,-0xA0(r1)                         PIPE
00A5CAB0 EB21FF60 ld r25,-0xA0(r1)                         70 (00A5CAAC) LHS[01]
00A5CAB4 D1A6000C stfs f13,0xC(r6)                          PIPE
00A5CAB8 D029000C stfs f1,0xC(r9)
00A5CABC D801FF60 stfd f0,-0xA0(r1)                         PIPE
00A5CAC0 E801FF60 ld r0,-0xA0(r1)                          70 (00A5CABC) LHS[01]
0

Pattern matching as an optimisation method, and a means of insight

2011-08-24T08:16:00.000+01:00

One feature of Haskell that I particularly like is its pattern-matching scheme. To route code down the correct codepath, Haskell uses pattern matching against data types and values, instead of extensive flow-control statements. This makes the code much simpler to read, as your code merely concerns itself with calling the correct functions, with the correct handler cases, and the underlying Haskell compiler deals with the flow-control minutiae.

I've increasingly found myself using a similar approach in C++ and HLSL using features such as polymorphic functions. A good example comes from optimisation: frequently, you find cases in the code where someone has written a "kitchen sink" function. This function internally handles many cases, and accordingly takes many parameters.

Frequently, many of these parameters are called with specificliteral values or constants. For example, many parameters may be zero which disable entire codepaths inside the function. For me, this is a bad code smell. But why? Many people would argue that since the code is branched away, what's the big deal? Especially if it means we can have just one definition of a function. Well, branches are very expensive on modern processors. On contemporary consoles, a floating point branch can yield a 27-cycle penalty. Not to mention all the preamble and branch set-up code.

I've therefore found myself increasingly providing multiple specialised versions of these functions, which omit some of the parameters that are permanently called with 0.0f or 1.0f and tailoring the code accordingly. Obviously this saves a lot of needless computation, branching and I$ misses.

What's interesting about this approach is that once you've performed a few of these refactorings, you start to expose patterns in the code and its structure that you previously were unaware of. You may find that your major axis of variation lies not in your chosen object hierarchy composition, but elsewhere. You may have a multitude of objects that add little value at a higher level, but now a large set of lower level functions that process the objects differently. Sometimes it makes more sense to structure and divide your code along the lines of what processing it does, not what you think the object "is", as that is the greater source of variation.

On the limitations of caching

2011-08-23T21:29:00.000+01:00

I've been spending a little time thinking about caching lately. Irradiance caching, optimising some caching in production code, and hardware caches. Here's a little brain-dump.

I've come to the conclusion that if you ever need anything more than a trivially small cache, your algorithm is wrong. This is particularly so for parallel programming, as caches are often a difficult to handle piece of interdependent common state.

My reasoning goes a little something like this. If you are relying on a cache for performance, it is usually to re-use the results of some expensive transformation a -> b. Crucially, there needs to be 1 a to one or more b (otherwise there is no cache reuse!). For an efficient cache, you need:

To be caching the results of a fairly complex computation,
To have an efficient means of storing the results of the computation,
To have a well-behaved cache indexing function.

Now, what often ends up happening is that you iterate over the many objects 'b', and populate the cache on-demand from a common ancestor piece of data, a. And here's the first problem. You are primarily iterating over the wrong dataset, and filtering and culling parts of the other dataset by proxy. You have an implicit data filtering or reduction operation that could be better handled more explicitly.

Ok, so, we switch our iteration round to iterate over the source objects 'a'. We iterate over each 'a' object, and produce the resulting 1 or more 'b' objects. At this point, do we really have a cache? We've turned our cache into a data transformation kernel. All of sudden our cache is no longer some intermediate state used in a computation, but it is now simply a piece of input data. Since we now have a piece of code transforming [a] -> [b], we have a piece of code that is very amenable to data parallelism. Since we've decoupled the code from a cache, it's also much easier to parallelise.

I think ultimately an extensive caching scheme is evidence that the data domain hasn't been adequately structured and managed. If you need a cache, you probably need more data conditioning, filtering, and reduction.

A good counter-example is the vertex stream processing unit of a modern GPU. These units read huge streams of data, transform them, and maintain minimal cached state. Vertex caches are typically miniscule in comparison to the volume of data they process. Yet a post-transform cache can often maintain a very high level of efficiency. The key to this cache's success is the pre-processing (ie structuring) of the data fed into it. It is not random. It is well conditioned, which yields efficiency at a number of levels, and enables a vertex cache to be used only as a relatively modest addition to the pipeline. The cache is the icing on the cake, not a fundamental component of the algorithm.

I think this is an important lesson when transitioning from writing single-threaded to parallel code. Don't dream up increasingly elaborate caching schemes. Instead think up separable, well-structured data-parallel transformations. When you feel the need for a cache, you probably need to carefully re-examine your data flow structure.

Where'd all the Haskell go?

2011-07-05T09:18:00.000+01:00

It's not gone. It'll be back.

I've recently moved house, and obviously that consumes a lot of time both before, and after the move. In order to write a Haskell article, I need to develop supporting code. I don't have much coding time right now.

On the other hand, it's pretty easy to knock out an article based on the knowledge gained during my Master's, and through my industry experience, so I'm writing on those subjects for a while.

Haskell will be back!

Photon Mapping in Haskell

2011-06-02T22:02:00.000+01:00

For a while I have been working on extending my Haskell raytracer to include a global illumination model. I chose to implement photon mapping, as it seemed to offer both a number of challenges, and it quite suited a functional programming approach.

Photon Tracing

One of the main challenges of implementing an algorithm like photon mapping is that it relies on a small amount of mutable state due to its monte-carlo integration. As each photon is traced through the scene, a random number is required to determine the fate of a photon for a given photon-to-material interaction. A random number generator typically requires state to iterate from one number to the next. Such state is easy to incorporate into imperative code, but it presents a challenge to the Haskell programmer.

There are two main challenges with state: management and propagation of state, and selecting what code should share a given state context. For example, in the photon mapper I could have introduced a single, global state for the random number generator. Or, I could have given each light source its own state. Or each photon. The greater scope that the state has, the greater the coupling in the code and the greater the difficulty in parallelisation.

I chose to give each photon traced its own random number generator state - and each photon had a different initial seed.

-- Realistic Image Synthesis Using Photon Mapping p60
tracePhoton :: [Photon] -> Photon -> SceneGraph -> StdGen -> (Int, Int) -> [Photon]
tracePhoton !currentPhotons (Photon !photonPower !photonPosDir) sceneGraph !rndState !(bounce, maxBounces) =
-- See if the photon intersects any surfaces
case findNearestIntersection sceneGraph ray of
Nothing -> currentPhotons
Just (obj, t, subId) -> case photonFate of
-- Diffuse reflection. Here, we store the photon that got reflected, and trace a new photon - but only if it's bright enough to be worthwhile
DiffuseReflect -> if Colour.magnitude newPhotonPower > brightnessEpsilon && (bounce + 1) <= maxBounces
then tracePhoton (storedPhoton : currentPhotons) reflectedPhoton sceneGraph rndState'' (bounce + 1, maxBounces)
else storedPhoton : currentPhotons
where
!reflectedPhoton = Photon newPhotonPower (surfacePos, reflectedDir)
!(reflectedDir, rndState'') = diffuseReflectionDirection rndState' tanSpace

-- Specular reflection. Here, we reflect the photon in the fashion that the surface would reflect towards the viewer and
-- aim to absorb it somewhere else in the photon map
SpecularReflect -> if Colour.magnitude newPhotonPower > brightnessEpsilon && (bounce + 1) <= maxBounces
then tracePhoton currentPhotons reflectedPhoton sceneGraph rndState' (bounce + 1, maxBounces)
else currentPhotons
where
!reflectedPhoton = Photon newPhotonPower (surfacePos, reflectedDir)
!reflectedDir = Vector.negate (snd photonPosDir) `reflect` normal

-- Absorb. The photon simply gets absorbed into the map
Absorb -> storedPhoton : currentPhotons
where
!(photonFate, rndState') = runState (choosePhotonFate coefficients) rndState
!coefficients = russianRouletteCoefficients (material obj)
!newPhotonPower = computeNewPhotonPower photonFate coefficients photonPower (material obj)
!tanSpace = primitiveTangentSpace (primitive obj) subId hitPosition obj
!normal = thr tanSpace
!hitPosition = pointAlongRay ray t
!surfacePos = hitPosition + (normal Vector.<*> surfaceEpsilon)
!brightnessEpsilon = 0.1
!storedPhoton = Photon photonPower (surfacePos, snd photonPosDir)
where
!ray = rayWithPosDir photonPosDir 10000

KD-Tree Construction

The resulting list of photons is inserted into a KD tree in order to efficiently locate the set of photons near to a point in space. This is a simple recursive operation:

buildKDTree :: [Photon] -> PhotonMapTree
buildKDTree photons
| length photons == 1 = PhotonMapLeaf (head photons)
| length photons >= 2 = let (boxMin, boxMax) = photonsBoundingBox photons
!axis = largestAxis (boxMax - boxMin)
!photonsMedian = foldr ((+) . fst . posDir) zeroVector photons Vector. fromIntegral (length photons)
!value = component photonsMedian axis
photonsGT = Prelude.filter (\x -> component ((fst . posDir) x) axis > value) photons
photonsLE = Prelude.filter (\x -> component ((fst . posDir) x) axis <= value) photons
in if length photonsGT > 0 && length photonsLE > 0
then PhotonMapNode axis value (buildKDTree photonsGT) (buildKDTree photonsLE)
else let (photons0', photons1') = trace "Using degenerate case" $ degenerateSplitList photons in PhotonMapNode axis value (buildKDTree photons0') (buildKDTree photons1')
| otherwise = error ("Invalid case, length of array is " ++ show (length photons) ++ "\n")

KD-Tree Query

The photon mapping algorithm uses an interesting combination of a kd tree and max-heap to locate only the closest N photons to a point. The max-heap is sorted on the squared distance to a photon. This means the most distant photons are stored towards the top of the max-heap, making them easy to discard.

Traversal of the tree commences with a point of interest and a radius. Any photons found within that radius are added to the max-heap. If the max-heap exceeds its specified size then excess, distant photons are dropped from the top of the heap. Since the max-heap stores the most distant photon at the top of the heap, it is easy to monitor the current-furthest photon and tighten our search radius during traversal.

-- Use a max heap to make it easy to eliminate distant photons
data GatheredPhoton = GatheredPhoton Float Photon deriving (Show)
type PhotonHeap = MaxHeap GatheredPhoton

instance Ord GatheredPhoton where
compare (GatheredPhoton dist1 _) (GatheredPhoton dist2 _)
| dist1 == dist2 = EQ
| dist1 <= dist2 = LT
| otherwise = GT

instance Eq GatheredPhoton where
(GatheredPhoton dist1 _) == (GatheredPhoton dist2 _) = dist1 == dist2

minimalSearchRadius !rSq !photonHeap = case viewHead photonHeap of
Nothing -> rSq
Just (GatheredPhoton dSq _) -> Prelude.min rSq dSq

-- Gather photons for irradiance computations
-- Algorithm adapted from Realistic Image Synthesis Using Photon Mapping p73
gatherPhotons :: PhotonMapTree -> Position -> Float -> PhotonHeap -> Int -> PhotonHeap
gatherPhotons (PhotonMapNode !axis !value gtChild leChild) !pos !rSq !photonHeap !maxPhotons
-- In this case, the split plane bisects the search sphere - search both halves of tree
| (value - posComponent) ** 2 <= rSq = let heap1 = gatherPhotons gtChild pos rSq' photonHeap maxPhotons
rSq'' = minimalSearchRadius rSq' heap1
heap2 = gatherPhotons leChild pos rSq'' photonHeap maxPhotons
newHeap = union heap1 heap2
in Data.Heap.drop (size newHeap - maxPhotons) newHeap

-- One side of the tree...
| posComponent > value = gatherPhotons gtChild pos rSq' photonHeap maxPhotons

-- ... or the other
| posComponent <= value = gatherPhotons leChild pos rSq' photonHeap maxPhotons

-- Prolapse
| otherwise = error "gatherPhotons: unexplained/unexpected case here"
where
!posComponent = component pos axis
!rSq' = minimalSearchRadius rSq photonHeap -- Refine search radius as we go down tree to search no further than closest allowed photon
gatherPhotons (PhotonMapLeaf !p) !pos !rSq !photonHeap !maxPhotons
| distSq < rSq = let newHeap = insert (GatheredPhoton distSq p) photonHeap
in Data.Heap.drop (size newHeap - maxPhotons) newHeap -- Discard any excess photons - we get rid of the furthest ones
| otherwise = photonHeap
where !distSq = pos `distanceSq` (fst . posDir) p

Parallelisation

The photon tracing stage is trivially easy to parallelise, provided that you have carefully treated the shared state required by the random number generator. Since I have provided each traced photon with it's own state, which is uniquely initialised, all photons are data independent and therefore easy to trace in parallel.

I have not yet parallelised the photon gathering stage, due to some issues with the current implementation's efficiency. Two options present themselves for parallelisation. The first is to parallelise the traversal of high-level sub-branches of the KD tree. This naive approach would yield some speedup, but the work of each thread is likely to be highly unbalanced.

Speculatively, a work-stealing queue type approach may help here. The high-level nodes to be traversed could be initially inserted into a queue. When a worker thread pulls data, it could compare how many nodes remain to be traversed against the total number of workers. If the ratio of nodes:workers becomes very low, nodes could be repeatedly removed from the queue and replaced with their child nodes, until the ratio becomes favourable. This would populate the queue with a fairly heterogenous mixture of large and small amounts of work to be done, which would help fill up and balance processor time across the threads.

Results

Here is the current output image, with all options cranked up:

The code requires refinement and tuning to eliminate the bias, improve the colour bleeding and fix a few minor artefacts.

The photon mapper is currently very slow. It can take hours on my MacBook Pro to trace a scene containing 100,000 photons, and that is a very small number of photons.

Most of the time is spent in the gathering phase. The current implementation is very space intensive and requires further serious attention. The photon emission and tracing step, however, is extremely quick, requiring only a few seconds to trace 100,000 photons.

Code

Full source code is available on github:

https://github.com/TomHammersley/HaskellRenderer

Future Work

I am now working to improve the efficiency of the gathering phase which will improve efficiency and development speed. I am also implementing an irradiance cache to further optimise execution speed. Clearly various minor bugs also need addressing.

References

Further implementation details can be found primarily in Henrik Wann Jensen's photon mapping book, supplemented by various SIGGRAPH course materials.

Haskell Renderer Source Code on Github

2011-05-15T07:16:00.001+01:00

Just a quick note to let anyone interested that the source code to my work-in-progress Haskell raytracer and photon mapper is now available at:

https://github.com/TomHammersley/HaskellRenderer

I'm currently mid way through debugging and optimising the photon map solution. It works, though it takes 10 minutes or so on my MacBook Pro.

When run it will output the following image:

To build the program you just need the latest Haskell Platform and the bmp package from cabal.

I'm still a beginner at Haskell so some of the code may not represent the most idiomatic or best practice Haskell. If anyone has any advice, comments or constructive criticism, I'm keen to hear it at tomhammersley@gmail.com.

(I plan to make the package Cabal compatible soon).

Cheers.

Designs That Leak State: Using Haskell Notation to Analyse C++ Design

2011-04-11T21:07:00.002+01:00

I came across an interesting design failure today. In short, a C++ object was internally storing the state specific to a piece of multi-threaded client code. The client code was leaking its state data into the C++ object. In generalised, abstracted terms:

class Object
{
Real member variables
...
unsigned char m_state1[MaxContexts];

float m_state2[MaxContexts];

bool m_state3[MaxContexts];

};

Each of the threads would read from 'Object', process it, and store the results in one of the elements of m_state.

It's a bad pattern, but I daresay a common pattern.

What is wrong here?

The first thing to notice is the series of member variables, all arrays, sized by MaxContexts. A series of variables like this is a dead giveaway that really there is semantically some concept in the system that is not explicitly represented in the code. Given that we are trying to represent an object in a series of "contexts", with different values for each, it suggests that these variables are not properties of the object, but actually variables of some other, auxiliary state. If there is logically one of something, it's a member variable of the class in question. If an object embeds a number of them, it's probably someone else's variable.

Fundamentally, my problem with this is that we're embedding some client's data within the object itself. When the client code operates on the object, it is no longer a pure transformation of one object into another, it is now a conceptually more complex, confusing transformation of one object, to the same object with new state.

This additional state presents serious problems to multi-threading a piece of code. When foreign state is embedded unnecesarily in an object like this it makes it much more difficult to manage the ownership, lifetime and correctness of data. All the other clients of the Object are burdened with helping manage the state of some other system. Appropriate ownership, coupling, and lifetime of data is all-important in parallel programming.

Our problem is that we have inverted the "has-a" relationship. The object does not "have a" state in a number of contexts. In reality, the contexts have the state and they "have", or refer to the object.

Haskell Function Signature Analysis

I increasingly find that Haskell's function signature notation is a lucid yet lightweight means of analysing data flow. There are many other alternative, successful methods, yet I find the Haskell analysis direct and natural. These days I find myself thinking directly in these terms to understand data transformations.

In Haskell terms, we have:

f :: Object -> Object

Or possible even:

f :: Object -> (Object, State)

We're not sure! We would prefer:

f :: Object -> State

Usage of State

The next question is, how is this state really used? It is not really part of an object's update function, as the multiple contexts suggest the data is not part of the object. So, what is the client code of the object doing with the state? It doesn't really matter. Whatever it's doing, we can infer that:
a) The state is a product of the object. If it were not, it would not be in the object.
b) We have multiple contexts, implying multiple passes over the data.

When a piece of state is the product of an object, it is likely that it is intermediate data. Given that we expect multiple passes over the data, we're going to be caching a lot of pieces of state data that have nothing to do with the context in question. We could switch the state representation from SoA to AoS, but it's still data in the wrong place.

Let's say we have 10 'contexts' and 100 objects. Let's say the object is 20 bytes, and the 'state' is 6 byte. We have a total data volume of 100 * (20 + 6 * 10) bytes, of which 100 * (20 + 6) bytes is used. Out of the total data volume, we're therefore using only 33% of the total data volume.

So, if we switch from the:

[Object] -> [Object]

implementation to a

[Object] -> [State]

implementation, we now have 100 * 20 + 100 * 6 bytes total data volume and 100% of the data is used.

Where is the data going?

The previous transformation is essentially a map operation. We're applying a function a list of Objects, and returning a list of State objects. A state is not an endpoint of a data processing; it implies that further processing is going to use the data.

Returning to Haskell, let's say our original state-producing function is f:

state = map f objects

Something else, g, is going to use that data:

output = map g (map f objects)

Or in terms of function application:

g . f

Here's our first big clue. All that data that originally lived inside the object, most likely on the heap, is really just intermediate data polluting the cache. It's just input to a succeeding function. We can most likely eliminate that "state" data altogether and pass the output directly through to the succeeding function. Or, we can store the state in some piece of intermediate thread-local storage. Either way, we've now got a clearer understanding of the real data flow in the program.

Eliminating the containing Object

The next question becomes, which bits of the Object really contribute to the new output State?

As the State is the product of an Object, and Object is non-trivial it is unlikely that the State is a product of the entire Object. It is probably some subset of member variables.

To use something of a straw man example, suppose the code in question is a visibility operation. Suppose that Object holds a bounding volume, some geometry, and the State holds various visibility test results, distances, LODs and the likely. Our transformation is therefore more likely to be:

BoundingVolume -> LodDescription -> State

Most of the Object's data is not used here. We've passed a subset of the Object's data to some function. This is a naive pattern in C++ code. Useful, general, common patterns of code are often hidden behind a facade of encapsulation. Unrelated pieces of data (in a given context) are brought in unnecessarily. When we develop a semi-abstract purely functional description of the operation, we can build a much clearer understanding of the data flow and dependencies.

Now imagine rewriting and structuring this operation in C++ based on this understanding of the data. Arguably an experienced programmer can arrive directly at the same conclusion, but I find the Haskell notation a useful intermediate tool for reasoning about data flows.

Ray-Triangle Intersection in Haskell

2011-04-07T08:21:00.000+01:00

Recently, I've been experimenting with different ray-triangle intersection algorithms. There are many alternatives out there, and many more optimisations and special cases, but I was looking for an version that ran quickly, and "suited" functional programming.

I settled on the half-planes solution. This algorithm defines a plane at each edge of the triangle by crossing the triangle's normal and the edge direction. To check if a triangle is inside the triangle, you simply check to see if a point (on the plane of the triangle) is on the same side of all of those planes.

This costs a little extra storage space, but it's a very fast test.

So. How does the code pan out? Well, pointInsideTriangle is simple enough:

distanceToPlane :: Shape -> Vector -> Float
distanceToPlane (Plane !norm !dist) !pos = pos `dot3` norm + dist
distanceToPlane _ _ = undefined

pointInsideTriangle :: Triangle -> Position -> Bool
pointInsideTriangle !tri !point = foldr (&&) True $ map (\pln -> (distanceToPlane pln point) >= 0) (halfPlanes tri)

This is best read right-to-left. The map applies a lambda function to each of the half planes. This function simply returns a bool, telling us if a point is on the correct side of a triangle. This yields a list of bools. We need to turn that into a single result; this is done using a fold to reduce many values to a single Bool. I've deliberately used a right-fold so that the fold lazily consumes only the necessary parts of the list. It can early-out.

So there we have it. One line of Haskell (with a couple of helper functions) to do the test.

The rest of the Haskell code for a ray-triangle intersection is pretty straightforward:

intersectRayTriangle :: Ray -> Object -> Triangle -> Bool -> Maybe (Float, Triangle)

intersectRayTriangle !ray !obj !triangle !doubleSided

| doubleSided == False && (direction ray) `dot3` (normal $ plane triangle) > 0 = Nothing

| otherwise = case shapeClosestIntersect (plane triangle) ray obj of

Nothing -> Nothing

Just (dist', _) -> if pointInsideTriangle triangle (pointAlongRay ray dist')

then Just (dist', triangle)

else Nothing

Here we have a simple guard condition to permit backface culling. If the triangle is not double-sided, we do a backface culling test and possibly reject it. If not, we fall through to the default case and intersect the ray against the plane, and then test that intersection point for containment in the triangle.

Intersecting against a list is an interesting case. A simple version would just map against the list, and then attempt to find the closest intersection out of the resulting list. However, if we use tail recursion, we can add a parameter to hold our current "state" and therefore maintain a current-closest intersection. This eliminates a second search-the-list step, and also permits an early-reject optimisation:

intersectRayTriangleList :: [Triangle] -> Int -> Maybe (Float, Int) -> Ray -> Object -> Maybe (Float, Int)

intersectRayTriangleList !(x:xs) !index !currentResult !currentRay !obj = intersectRayTriangleList xs (index + 1) newResult newRay obj

where

(newRay, newResult) = case intersectRayTriangle currentRay obj x False of

Nothing -> (currentRay, currentResult)

Just (dist, _) -> (shortenRay currentRay dist, Just (dist, index))

intersectRayTriangleList [] _ !currentResult _ _ = currentResult

The pattern of "searching a list of x and return the closest" occurs frequently in raytracing. This is a potential candidate to be factored out into a re-usable function, or possibly even a monad.

The resulting code is also quite efficient. I'm now running at less than a minute for my test scene of the Cornell Box with 8x8 distributed raytracing. This is 30% quicker than my previously Möller-Trumbore ray-triangle intersection test.

How to Parallelise A Haskell Raytracer

2011-04-03T22:28:00.003+01:00

One of the big attractions of functional languages is that they're theoretically very easily to parallelise as they have no problematic side-effects You only describe how to perform a calculation, and never touch on the details of moving data around the machine. The run-time is then free to process data, given your functional expressions, in parallel.

I've been experimenting with Haskell's basic parallelism primitives, par, and pseq. These primitives trigger ("spark") evaluation of expressions in parallel - but only if the runtime deems it necessary.

After many failed or mildly successful experiments, I've finally managed to get a speedup. In short, I stopped trying to reinvent the wheel myself. I simply used parListChunk:

rayTraceImage :: Camera -> Int -> Int -> SceneGraph -> [Light] -> [Colour]
rayTraceImage camera renderWidth renderHeight sceneGraph lights = map (clamp . tracePixel eyePosition sceneGraph lights) rayDirections `using` parListChunk 256 rseq
where !rayDirections = [makeRayDirection renderWidth renderHeight camera (x, y) | y <- [0..(renderHeight - 1)], x <- [0..(renderWidth - 1)]]
!eyePosition = Camera.position camera

Here, we break the list of rays to be traced into data parallel chunks of 256, each free to be evaluated by a different thread. Each of these chunks is processed using map. Map applies the raytracing function to convert a ray direction into a colour, and we clamp the result.

And that's it.

Running it on two threads on my dual-core MacBook Pro halves the execution time.

What's nice is that the core 'map' expression is separated from the details of parallelisation. We've separated the concepts of raytracing and its parallel decomposition into two separate concerns. We are therefore free to vary either without intertwining their concerns.

Not a critical section, lock-free queue or atomic operation in sight!

Clearly, the nature of the functional approach surrenders a lot of low-level implementation details to your platform's runtime. However, it is more frequently the case that when you require such low-level intervention, you are working around deficiencies in your run-time platform, rather than dealing with fundamental algorithm details. Parallelising many algorithms follows established patterns.

It is also worth noting that the raytracer code is currently pretty straightforward code. I've not implemented any particular tricky optimisations, and certainly none of the ideas in work such as Ingo Wald's parallel raytracing PhD thesis. Not yet, anyway...

Laziness vs Strictness in Haskell

2011-03-24T13:38:00.001+00:00

I've just taken over a minute of execution time off my raytracer.

One of Haskell's biggest language features is the lazy evaluation system. By default, Haskell evaluates nearly all expressions lazily, as needed. This is one of the features that drew me to Haskell, and indeed it seems a reasonably apt way of writing a raytracer. Raytracing often consists of selectively traversing different parts of large databases for different parts of the scene. Many acceleration schemes involve simple mechanisms to lazily evaluate results.

Laziness, however, does pose an overhead. Haskell introduces a "thunk" for every lazy expression. A thunk describes how to calculate a deferred expression when required. For some expressions, a thunk is a potential time saver: you can defer work until needed, calculate the required value, then re-use it. For other things, such as simple arithmetic operations like vector addition, a thunk is an overhead often comparable to the original operation.

Haskell offers a couple of ways to avoid thunks and therefore lazy evalation. The main two I use are the BangPatterns extension to explicitly mark something as strictly (ie, non-lazily) evaluated, and seq.

The BangPatterns extension allows you to mark expressions that you want to be explicitly evaluated with a !. So, for a Vector, you would write Vector !Float !Float !Float !Float. This way, you force anything you assign to a vector to be strictly evaluated and stored in the vector.

Seq is used to strictly force the evaluation of a function. Seq has the signature of a -> b -> a. It therefore takes two parameters, and returns the second.

An interesting example of where lazy evaluation proven to be an overheard, rather than a win was my light accumulation function. Given a point in space to be shaded, this accumulates the contribution of every light source. Now, the inputs to this function are going to be different every time this function is called as every ray intersection is going to be unique. Thunks will therefore present an additional overhead.

I rewrote my function to use seq. I use seq to strictly evaluate the application of this function, then process the remainder of the list:

accumulateLight :: [Light] -> Colour -> SceneGraph -> Vector -> Vector -> Material -> Vector -> Colour
accumulateLight (x:xs) acc sceneGraph shadePos normal objMaterial viewDirection = let result = acc + (applyLight sceneGraph shadePos normal objMaterial viewDirection x)
in seq result (accumulateLight xs result sceneGraph shadePos normal objMaterial viewDirection)
accumulateLight [] acc _ _ _ _ _ = acc

This change alone saved over a minute of execution time! Since a raytracer generally consists of many similar kinds of loops, I am now searching for wider application of this technique. It's a very difficult call to make with a limited data set. Some functions may be faster if strictly evaluated with a small dataset, but may perform better if evaluated lazily with a larger dataset.

Profiling will guide me.

Pragmatism vs Idealism (DOD vs OOP)

2011-03-23T09:52:00.002+00:00

So, it seems no blog these days is complete without a yet another stab at the DOD vs OOP argument.

Here's my take.

I could have picked any number of dichotomy based titles for this post. Objective vs Subjective. Pragmatism vs Idealism. Design vs Performance. PC vs Console. There is a tendency to view DOD and OOP as somewhat polar extremes that you are either for, or against. Personally, I see merit in both camps. It is something of a false dichotomy. I see value in the ideals of OOP and DOD. Equally, I see shortcomings in both camps. For me, it's not an either/or, it's a question of how you blend them together.

For me, the real key issue is choosing what approach to adopt for a given task. This is where the real challenge lies.

The arguments used to support OOP are essentially based on intangible qualities. Design. Re-use. Extensibility. Maintainability. Elegance. Whilst valuable, these are generally somewhat subjective qualities. What constitutes good design is based on widely held best practices, and those practices are based on a consensus of constantly evolving expert opinion. It's very hard to measure these things explicitly. It's very hard to say how much "better" a more OOP centric approach can make a piece of code.

The arguments for DOD, by contrast, are based on real-world, objective data: performance. You can directly measure performance, and directly state the improvement that a DOD approach has yielded for a given task. Similarly, you can measure the cost of an OOP based approach in terms of lost performance. Performance has a useful, marginal benefit in a software project. Other qualities have their benefits but are frequently less directly felt.

You therefore end up with a challenging situation. How can the intangible qualities that support OOP argue against the real-world, objective data supporting DOD? How can we support OOP over DOD when it is based on subjective preferences, leading to design decisions that are somewhat arbitrary. Especially when you are trying to a ship a console game!

As I see it, the performance data that supports DOD should be supporting both arguments. The current generation of consoles in many ways are indicative of what future computing platforms will generally look like. Many cores. Slow memory. Long, simple instruction pipelines.

OOP has its positive qualities, but it is important to view it in context: it largely emerged as a means of improving the structure of large software systems by breaking software into singular, encapsulated objects. Whilst these concerns are significant, there are other, equally valid ways of addressing those concerns than OOP alone. And with the continuing, accelerating trend of parallel hardware, OOP alone cannot be positioned as the primary concern in designing software systems.

For me, my primary concern in software design is parallelism first, DOD second and OOP third.

Sphere Tree Construction and Traversal in Haskell

2011-03-22T19:33:00.001+00:00

Finally. I've finally finished writing and debugging my sphere tree construction code in Haskell.

It should have been so easy. Tree structures are very easily recursively defined which is ideal for functional languages like Haskell. Take your data, split it up, and build the children by recursion. Simple. Only it was anything but...

The logic of the task was simple enough. The main problem was debugging. This is the sort of task that is prone to simple mistakes that cause time-consuming debugging problems. When working in an imperative language, you can easily put down a breakpoint and inspect your data, or if you recur too much the stack blows up, or you can step through your code and see why you're getting stuck in loops. These strategies don't map well to Haskell.

What does work in Haskell is to start small, build up gradually and constantly verify your results. You end up with a development cycle that is fundamentally very similar to test driven development. You want to prove your code one little function at a time by feeding it controlled inputs and inspecting the outputs. In Haskell, your primary debugger is your brain - you need to think very carefully through your code and how your express your algorithms. Functional languages reduce the scope of creating simple bugs due to side effects, and Haskell in particular eliminates a lot of plumbing and housekeeping code. What's left is purely your algorithm - inspect it carefully, that is where your problem lies.

Bounding Volume Hierarchies

I considered a few approaches for my bounding volume hierarchy, and for now, I've settled on a fairly generic implementation of bounding volume hierarchies that currently is configured as a sphere tree.

Each tree node can optionally (Maybe) hold an object. Each node has a sphere, but that could later be generalised. The children are defined as a list, allowing me to switch using a quadtree, kdtree or octree type approach.

(Note: The code below is best copied and viewed in a real text editor. Also, I am not yet a Haskell expert, and there may be many superior ways to achieve these results. My purpose in writing is to share my experiences as a learner with other learners.)

data SphereTreeNode = SphereTreeNode { object :: Maybe Object, children :: [SphereTreeNode], boundingRadius :: Float, boundingCentre :: Vector } deriving (Show, Read)

Sphere Tree Construction

The tree itself is constructed using a fairly simple piece of Haskell:

-- Build up a sphere tree
buildSphereTree :: ([Object] -> [[Object]]) -> [Object] -> SphereTreeNode
buildSphereTree _ (obj : []) = SphereTreeNode (Just obj) [] nodeRadius nodeCentre
   where
   nodeCentre = calculateMeanPosition (obj:[])
   nodeRadius = calculateBoundingRadius (obj:[]) nodeCentre
buildSphereTree builder (obj:objs)
   | length (obj:objs) == 1 = error "Should have been handled by a different pattern"
   | length (obj:objs) == 0 = error "Should not have zero objects"
   | otherwise = SphereTreeNode Nothing nodeChildren nodeRadius nodeCentre
   where
   nodeCentre = calculateMeanPosition (obj:objs)
   nodeRadius = calculateBoundingRadius (obj:objs) nodeCentre
   nodeChildren = map (buildSphereTree builder) (builder (obj:objs))
buildSphereTree _ [] = error "Should not hit this pattern for buildSphereTree"

-- Find the mean of a collection of objects

calculateMeanPosition' :: [Object] -> Vector -> Vector

calculateMeanPosition' (obj : objects) acc = calculateMeanPosition' objects acc + (getCentre obj)

calculateMeanPosition' [] acc = acc

calculateMeanPosition :: [Object] -> Vector

calculateMeanPosition objects = setWTo1 ((calculateMeanPosition' objects zeroVector) fromIntegral (length objects))

-- Find the overall bounding radius of a list of objects

calculateBoundingRadius :: [Object] -> Vector -> Float

calculateBoundingRadius objs centre = foldr Prelude.max 0 (map (\obj -> shapeBoundingRadius (shape obj) (transform obj) centre) objs)

This object simply recursively builds a tree until it encounters an object list of size 1.

I've deliberately exercised a couple of Haskell idioms here. I've used another tail-recursive loop to calculate the mean centre of the objects, and a map/fold pair to calculate the sphere's bounding radius.

The most interesting part is the first parameter to buildSphereTree. I pass in a function of type ([Object] -> [[Object]]). This function is responsible for dividing the object list into a number of lists, each of which is the object list used to build a new child node. This trick allows me to abstract out the specific algorithm for building the tree into a user-supplied function.

In my current code, I'm using a simple KD-type approach:

-- Generate a plane to split the objects along

makeSplittingPlane :: AABB -> (Vector, Float)

makeSplittingPlane (boxMin, boxMax) = case largestAxis (boxMax - boxMin) of

                                        0 -> (xaxis, -(vecX midPoint))

                                        1 -> (yaxis, -(vecY midPoint))

                                        2 -> (zaxis, -(vecZ midPoint))

                                        _ -> error "Undefined value"

    where

      midPoint = (boxMin + boxMax) <*> 0.5

-- Make children using a kd tree

generateSceneGraphUsingKDTree :: [Object] -> [[Object]]

generateSceneGraphUsingKDTree objs = leftObjects : rightObjects : []

    where

      objBox = objectListBoundingBox objs

      (planeNormal, planeDist) = makeSplittingPlane objBox

      (leftObjects, rightObjects) = partition (\obj -> planeDist + (dot3 planeNormal $ getCentre obj) > 0) objs

We simply make a splitting plane and split the list of objects according to what side of the plane their centre lies: [[Object], [Object]].

I could easily replace this function with an alternative strategy. Or even, make a function that returns a function based on the characteristics of the presented list...

Bounding Volume Hierarchy Intersection

Now this is the real meat of the implementation. The function again takes the form of a tail-recursive loop. The code is essentially similar to implementing a recursive function using a software stack in C. The tail recursion is used to make a simple loop-type recursion rather than a true state-accumulation call-type recursion.

The code is somewhat complex as several optimisations have been added.

The code is passed a list of nodes to process. We process the first node in the list, and recur with the remainder of the list - and possibly, some newly added nodes. The code seeks to find and maintain the closest-intersecting object.

The code constructs a tuple of three values - a result for this object, a list of additional nodes to process and a ray. If an object or bounding volume is hit, this tuples holds new values. If no intersection is found, we simply continue with the existing state passed into the function.

The code first considers whether the given node has any children; if it has none, it only makes sense to test the object that the node should hold. We report an error if the node has no object (a node without children or objects is useless!).

If the node does have children, it holds a subtree we wish to accept or reject. We first test the bounding volume. If this is rejected, we simply retain the current results of the loop with no additional nodes to process. If the bounding volume is intersected, then we will at least recur to the children nodes. If the node contains an object and it is intersected, we update our current intersection results and shorten our ray.

Finally, we pattern match the empty array case to simply return the current results.

Here is the code:

intersectSphereTree :: [SphereTreeNode] -> Ray -> Maybe (Object, Float, Int) -> Maybe (Object, Float, Int)

intersectSphereTree (node:nodes) ray currentHit = intersectSphereTree (newNodeList ++ nodes) newRay thisResult

where

-- Intersect the ray with the bounding volume of this node

(thisResult, newNodeList, newRay) = case children node of

-- If the node has no children, don't bother with it's bounding volume and just check the object (if it has one)

[] -> case object node of

Nothing -> error "A node with no children should hold an object"

Just obj -> case shapeIntersect (shape obj) ray (transform obj) of

-- Didn't hit the object. Retain the current hit, and continue with remaining nodes on the list

Nothing -> (currentHit, [], ray)

-- We did hit this object. Update the intersection, and continue with remaining nodes on the list

Just (objHitDistance, objHitId) -> (Just (obj, objHitDistance, objHitId), [], shortenRay ray objHitDistance)

-- We have children. In this case it makes sense to test our bounding volume

nodeChildren -> case shapeIntersect (Sphere (boundingRadius node)) ray (translationMatrix' (boundingCentre node)) of -- (make a sphere centred at the object's transform matrix with given radius)

-- If we do not find an intersection, we do not update the results and we offer no further nodes to be traversed, thus skipping this subtree

Nothing -> (currentHit, [], ray)

-- If we do find an intersection against the bounding volume, then we try again against the actual object (if present)

Just (_, _) -> case object node of

Nothing -> (currentHit, nodeChildren, ray) -- No object; just pass to the children

Just obj -> case shapeIntersect (shape obj) ray (transform obj) of

-- Didn't hit the object. Retain the current hit, but offer up the children of the node as we hit the bounding volume

Nothing -> (currentHit, nodeChildren, ray)

-- We did hit this object. Update the intersection, and continue with the bounding volume's children

Just (objHitDistance, objHitId) -> (Just (obj, objHitDistance, objHitId), nodeChildren, shortenRay ray objHitDistance)

intersectSphereTree [] _ currentHit = currentHit

The striking thing about this code is that it steadily marches to the right of the screen, with each indentation being a new test of success or failure. Ideally, this code could be rewritten as a monad.

What About Infinite Objects?

They're simply "filter"ed out at the time of the sphere tree construction, and added to a separate list. It makes little sense to include an infinite object into any bounding volume hierarchy.

Results

The code continues to render my test scene as expected. Disappointingly, it is currently a little slower; however, this is not due to the scene graph code. I have recently rewritten my vector classes to be instances of the Num typeclass, so that I can simply write a + b, rather than coining new operators such as a <+> b. I haven't yet fully optimised my new code.

Other Observations

As I get further into Haskell, I'm finding myself increasingly thinking in terms of Lambda calculus. I no longer find myself thinking of what iterative steps I need to do to complete an operation, rather I find myself thinking in terms of the transformations of the underlying data and accompanying functions.

I find this new perspective very helpful when returning to other languages. It gives you a somewhat stream-processing-like perspective and thought process when facing standard C++ code.

What's Next?

The main focus now is to diagnose, understand and resolve the issues with my new typeclass-based fundamental data types such as Vector and Colour. I want to get it below the ~70 seconds it was previously taking. Following that, it is now time to build a parser and load some standard test scenes that will allow me to develop and optimise the ray tracer further.

A worked Haskell example

2011-03-03T20:59:00.001+00:00

I thought it might be interesting to do a worked example of some work-in-progress Haskell to discuss what I like about it.

The case in question is the code that traverses the scene graph and intersects a ray.

Now, I'm putting the ultimate efficiency and DOD type concerns aside for a moment. I grant you there may be more efficient/better ways of doing this, it's just for discussion purposes.

Disclaimer: I am not (yet) a Haskell expert. There are probably better ways of doing this or writing it in Haskell. I write for the reasons of discussing what I've learnt about trying to write a practical application in Haskell. If there are better ways of doing it, I'd love to know, as I'm still learning.

So, first up, the code:

intersectSceneGraph' :: [SceneGraphNode] -> Ray -> Maybe (Object, Position, Int) -> Maybe (Object, Position, Int)
intersectSceneGraph' (node:nodes) ray currentHit = intersectSceneGraph' (newNodeList ++ nodes) ray thisResult
   where
   (thisResult, newNodeList) = case shapeIntersect (Sphere $ boundingRadius node) ray (transform (object node)) of
   Nothing -> (currentHit, [])
   Just (_, _) -> case shapeIntersect nodeObjectShape ray nodeObjectTransform of
   Nothing -> (currentHit, children node)
   Just (objHitDistance, objHitId) -> (Just (object node, pointAlongRay ray objHitDistance, objHitId), children node)
   where
   nodeObjectShape = shape $ object node
   nodeObjectTransform = transform $ object node
intersectSceneGraph' [] _ currentHit = currentHit

intersectSceneGraph :: SceneGraphNode -> Ray -> Maybe (Object, Position, Int)
intersectSceneGraph node ray = intersectSceneGraph' (node : []) ray Nothing

(If the code does not display clearer in your browser, cut and paste to a text editor)

This function is split into two parts. intersectSceneGraph is essentially a front end, a wrapper for intersectSceneGraph'. intersectSceneGraph' traverses a list of scene graph nodes, tests the ray against the bounding volume, and if it intersets it, it then intersects against the contained object and recurs to the children.

So, how does it traverse a tree?

The Haskell code is very similar to C code that traverses a tree using a software stack.

The operation is defined in a tail recursive fashion. I process one element of the list, and then recur to the remainder of the list. If the node in question is a new closest hit, this gets passed through to the next recursion as the new current state. If the bounding volume is intersected, then I add the children of the node to the list. If it is not intersected, I add nothing, an empty list.

This business of passing through state seems a little odd. Why would one instance of a function pass its internal state to another? This is a Haskell idiom. By passing the state from one call to another, as an "accumulator" parameter, it allows the compiler to eliminate any intermediate state on the stack and implement the recursion as a goto loop.

And all of this is lazily evaluated.

In pseudo-C, the code might look something like this:

closestHit = null;

while((topOfStack = stack.pop()) != NULL)
{

if(intersect(ray, topOfStack->bounding volume))

{

if(intersect(ray, topOfStack->object))

{

closestHit = this one

}

push(stack, topOfStack->children);

}
}

I think this shows the elegance of Haskell. Think of what is not there. There are no braces to explain scope to the compiler: you just indent it. There are no incidental variables describing the state of loop variables and so on. The compiler automatically pattern matches and dispatches to appropriate special cases instead. No pointer walks. No null checks. No explicit shuffling of data from A to B. You have a closer, more direct description of your algorithm, and the compiler fills in the gubbins.

Of course, there is a cost to all of this. By sacrificing that lower-level control, you lose out on the ultimate optimisability of the code. Still, the objective here is to find the best way to learn and use Haskell, rather than write the ultimate real-time raytracer in SIMD assembly! :-)