And what happened?

Fast 2D and 3D Hilbert curves and Morton codes

2011-08-31T19:07:00.003+01:00

Hilbert curves, Morton codes ad other space filling algorithms can be used to improve cache coherency, disc access and general sorting and searching.

Presented here are some compact, non-recursive bit twiddling C implementations for fast integer conversions between coordinates and curve indices.

While the Morton code to Hilbert curve conversion functions do still have loops in them, these can be unrolled if the required bit counts are known.

The 'bits' parameter in the Morton/Hilbert conversion functions refers to the extent of the encoded space. For instance a bits value of 1 indicates a 2x2(x2) space, 2 is a 4x4(x4) space, 3 is a 8x8(x8) space and so on.

Note: uint is assumed to be a 32-bit unsigned integer.

    uint MortonToHilbert2D( const uint morton, const uint bits )
    {
      uint index = 0;
      uint remap = 0xb4;
      uint block = ( bits << 1 );
      while( block )
      {
        block -= 2;
        uint mcode = ( ( morton >> block ) & 3 );
        uint hcode = ( ( remap >> ( mcode << 1 ) ) & 3 );
        remap ^= ( 0x82000028 >> ( hcode << 3 ) );
        index = ( ( index << 2 ) + hcode );
      }
      return( index );
    }
    
    uint HilbertToMorton2D( const uint hilbert, const uint bits )
    {
      uint index = 0;
      uint remap = 0xb4;
      uint block = ( bits << 1 );
      while( block )
      {
        block -= 2;
        uint hcode = ( ( hilbert >> block ) & 3 );
        uint mcode = ( ( remap >> ( hcode << 1 ) ) & 3 );
        remap ^= ( 0x330000cc >> ( hcode << 3 ) );
        index = ( ( index << 2 ) + mcode );
      }
      return( index );
    }
    
    uint MortonToHilbert3D( const uint morton, const uint bits )
    {
      uint index = morton;
      if( bits > 1 )
      {
        uint block = ( ( bits * 3 ) - 3 );
        uint hcode = ( ( index >> block ) & 7 );
        uint mcode, shift, signs;
        shift = signs = 0;
        while( block )
        {
          block -= 3;
          hcode <<= 2;
          mcode = ( ( 0x20212021 >> hcode ) & 3 );
          shift = ( ( 0x48 >> ( 7 - shift - mcode ) ) & 3 );
          signs = ( ( signs | ( signs << 3 ) ) >> mcode );
          signs = ( ( signs ^ ( 0x53560300 >> hcode ) ) & 7 );
          mcode = ( ( index >> block ) & 7 );
          hcode = mcode;
          hcode = ( ( ( hcode | ( hcode << 3 ) ) >> shift ) & 7 );
          hcode ^= signs;
          index ^= ( ( mcode ^ hcode ) << block );
        }
      }
      index ^= ( ( index >> 1 ) & 0x92492492 );
      index ^= ( ( index & 0x92492492 ) >> 1 );
      return( index );
    }
    
    uint HilbertToMorton3D( const uint hilbert, const uint bits )
    {
      uint index = hilbert;
      index ^= ( ( index & 0x92492492 ) >> 1 );
      index ^= ( ( index >> 1 ) & 0x92492492 );
      if( bits > 1 )
      {
        uint block = ( ( bits * 3 ) - 3 );
        uint hcode = ( ( index >> block ) & 7 );
        uint mcode, shift, signs;
        shift = signs = 0;
        while( block )
        {
          block -= 3;
          hcode <<= 2;
          mcode = ( ( 0x20212021 >> hcode ) & 3 );
          shift = ( ( 0x48 >> ( 4 - shift + mcode ) ) & 3 );
          signs = ( ( signs | ( signs << 3 ) ) >> mcode );
          signs = ( ( signs ^ ( 0x53560300 >> hcode ) ) & 7 );
          hcode = ( ( index >> block ) & 7 );
          mcode = hcode;
          mcode ^= signs;
          mcode = ( ( ( mcode | ( mcode << 3 ) ) >> shift ) & 7 );
          index ^= ( ( hcode ^ mcode ) << block );
        }
      }
      return( index );
    }

    uint Morton_2D_Encode_5bit( uint index1, uint index2 )
    { // pack 2 5-bit indices into a 10-bit Morton code
      index1 &= 0x0000001f;
      index2 &= 0x0000001f;
      index1 *= 0x01041041;
      index2 *= 0x01041041;
      index1 &= 0x10204081;
      index2 &= 0x10204081;
      index1 *= 0x00108421;
      index2 *= 0x00108421;
      index1 &= 0x15500000;
      index2 &= 0x15500000;
      return( ( index1 >> 20 ) | ( index2 >> 19 ) );
    }

    void Morton_2D_Decode_5bit( const uint morton, uint& index1, uint& index2 )
    { // unpack 2 5-bit indices from a 10-bit Morton code
      uint value1 = morton;
      uint value2 = ( value1 >> 1 );
      value1 &= 0x00000155;
      value2 &= 0x00000155;
      value1 |= ( value1 >> 1 );
      value2 |= ( value2 >> 1 );
      value1 &= 0x00000133;
      value2 &= 0x00000133;
      value1 |= ( value1 >> 2 );
      value2 |= ( value2 >> 2 );
      value1 &= 0x0000010f;
      value2 &= 0x0000010f;
      value1 |= ( value1 >> 4 );
      value2 |= ( value2 >> 4 );
      value1 &= 0x0000001f;
      value2 &= 0x0000001f;
      index1 = value1;
      index2 = value2;
    }

    uint Morton_2D_Encode_16bit( uint index1, uint index2 )
    { // pack 2 16-bit indices into a 32-bit Morton code
      index1 &= 0x0000ffff;
      index2 &= 0x0000ffff;
      index1 |= ( index1 << 8 );
      index2 |= ( index2 << 8 );
      index1 &= 0x00ff00ff;
      index2 &= 0x00ff00ff;
      index1 |= ( index1 << 4 );
      index2 |= ( index2 << 4 );
      index1 &= 0x0f0f0f0f;
      index2 &= 0x0f0f0f0f;
      index1 |= ( index1 << 2 );
      index2 |= ( index2 << 2 );
      index1 &= 0x33333333;
      index2 &= 0x33333333;
      index1 |= ( index1 << 1 );
      index2 |= ( index2 << 1 );
      index1 &= 0x55555555;
      index2 &= 0x55555555;
      return( index1 | ( index2 << 1 ) );
    }
    
    void Morton_2D_Decode_16bit( const uint morton, uint& index1, uint& index2 )
    { // unpack 2 16-bit indices from a 32-bit Morton code
      uint value1 = morton;
      uint value2 = ( value1 >> 1 );
      value1 &= 0x55555555;
      value2 &= 0x55555555;
      value1 |= ( value1 >> 1 );
      value2 |= ( value2 >> 1 );
      value1 &= 0x33333333;
      value2 &= 0x33333333;
      value1 |= ( value1 >> 2 );
      value2 |= ( value2 >> 2 );
      value1 &= 0x0f0f0f0f;
      value2 &= 0x0f0f0f0f;
      value1 |= ( value1 >> 4 );
      value2 |= ( value2 >> 4 );
      value1 &= 0x00ff00ff;
      value2 &= 0x00ff00ff;
      value1 |= ( value1 >> 8 );
      value2 |= ( value2 >> 8 );
      value1 &= 0x0000ffff;
      value2 &= 0x0000ffff;
      index1 = value1;
      index2 = value2;
    }

    uint Morton_3D_Encode_5bit( uint index1, uint index2, uint index3 )
    { // pack 3 5-bit indices into a 15-bit Morton code
      index1 &= 0x0000001f;
      index2 &= 0x0000001f;
      index3 &= 0x0000001f;
      index1 *= 0x01041041;
      index2 *= 0x01041041;
      index3 *= 0x01041041;
      index1 &= 0x10204081;
      index2 &= 0x10204081;
      index3 &= 0x10204081;
      index1 *= 0x00011111;
      index2 *= 0x00011111;
      index3 *= 0x00011111;
      index1 &= 0x12490000;
      index2 &= 0x12490000;
      index3 &= 0x12490000;
      return( ( index1 >> 16 ) | ( index2 >> 15 ) | ( index3 >> 14 ) );
    }

    void Morton_3D_Decode_5bit( const uint morton,
                                uint& index1, uint& index2, uint& index3 )
    { // unpack 3 5-bit indices from a 15-bit Morton code
      uint value1 = morton;
      uint value2 = ( value1 >> 1 );
      uint value3 = ( value1 >> 2 );
      value1 &= 0x00001249;
      value2 &= 0x00001249;
      value3 &= 0x00001249;
      value1 |= ( value1 >> 2 );
      value2 |= ( value2 >> 2 );
      value3 |= ( value3 >> 2 );
      value1 &= 0x000010c3;
      value2 &= 0x000010c3;
      value3 &= 0x000010c3;
      value1 |= ( value1 >> 4 );
      value2 |= ( value2 >> 4 );
      value3 |= ( value3 >> 4 );
      value1 &= 0x0000100f;
      value2 &= 0x0000100f;
      value3 &= 0x0000100f;
      value1 |= ( value1 >> 8 );
      value2 |= ( value2 >> 8 );
      value3 |= ( value3 >> 8 );
      value1 &= 0x0000001f;
      value2 &= 0x0000001f;
      value3 &= 0x0000001f;
      index1 = value1;
      index2 = value2;
      index3 = value3;
    }

    uint Morton_3D_Encode_10bit( uint index1, uint index2, uint index3 )
    { // pack 3 10-bit indices into a 30-bit Morton code
      index1 &= 0x000003ff;
      index2 &= 0x000003ff;
      index3 &= 0x000003ff;
      index1 |= ( index1 << 16 );
      index2 |= ( index2 << 16 );
      index3 |= ( index3 << 16 );
      index1 &= 0x030000ff;
      index2 &= 0x030000ff;
      index3 &= 0x030000ff;
      index1 |= ( index1 << 8 );
      index2 |= ( index2 << 8 );
      index3 |= ( index3 << 8 );
      index1 &= 0x0300f00f;
      index2 &= 0x0300f00f;
      index3 &= 0x0300f00f;
      index1 |= ( index1 << 4 );
      index2 |= ( index2 << 4 );
      index3 |= ( index3 << 4 );
      index1 &= 0x030c30c3;
      index2 &= 0x030c30c3;
      index3 &= 0x030c30c3;
      index1 |= ( index1 << 2 );
      index2 |= ( index2 << 2 );
      index3 |= ( index3 << 2 );
      index1 &= 0x09249249;
      index2 &= 0x09249249;
      index3 &= 0x09249249;
      return( index1 | ( index2 << 1 ) | ( index3 << 2 ) );
    }
    
    void Morton_3D_Decode_10bit( const uint morton,
                                 uint& index1, uint& index2, uint& index3 )
    { // unpack 3 10-bit indices from a 30-bit Morton code
      uint value1 = morton;
      uint value2 = ( value1 >> 1 );
      uint value3 = ( value1 >> 2 );
      value1 &= 0x09249249;
      value2 &= 0x09249249;
      value3 &= 0x09249249;
      value1 |= ( value1 >> 2 );
      value2 |= ( value2 >> 2 );
      value3 |= ( value3 >> 2 );
      value1 &= 0x030c30c3;
      value2 &= 0x030c30c3;
      value3 &= 0x030c30c3;
      value1 |= ( value1 >> 4 );
      value2 |= ( value2 >> 4 );
      value3 |= ( value3 >> 4 );
      value1 &= 0x0300f00f;
      value2 &= 0x0300f00f;
      value3 &= 0x0300f00f;
      value1 |= ( value1 >> 8 );
      value2 |= ( value2 >> 8 );
      value3 |= ( value3 >> 8 );
      value1 &= 0x030000ff;
      value2 &= 0x030000ff;
      value3 &= 0x030000ff;
      value1 |= ( value1 >> 16 );
      value2 |= ( value2 >> 16 );
      value3 |= ( value3 >> 16 );
      value1 &= 0x000003ff;
      value2 &= 0x000003ff;
      value3 &= 0x000003ff;
      index1 = value1;
      index2 = value2;
      index3 = value3;
    }

Row major and column major matrices

2011-05-22T17:40:00.004+01:00

Row major and column major matrices (or proving that black is white)

There's a good chance that you already think you understand the difference between row major and column major matrices, but are you right?

Let's go over the standard definition of row and column major:

A matrix is row major if the elements each row are contiguous in memory and column major if the elements of each column are contiguous in memory.

For example, when viewed as a 1 dimensional array, the elements of a 4 by 4 row major matrix are ordered:
r1c1, r1c2, r1c3, r1c4, r2c1, r2c2, r2c3, r2c4, r3c1, r3c2, r3c3, r3c4, r4c1, r4c2, r4c3, r4c4

And the elements a 4 by 4 column major matrix are ordered:
r1c1, r2c1, r3c1, r4c1, r1c2, r2c2, r3c2, r4c2, r1c3, r2c3, r3c3, r4c3, r1c4, r2c4, r3c4, r4c4

In a row major matrix the step between columns is 1 element (a small or minor step) and the step between rows is [column] elements (a larger or major step).

In a column major matrix the step between rows is 1 element (a small or minor step) and the step between columns is [row] elements (a larger or major step).

Viewed as a 1 dimensional array, a column major matrix is the same as the result of applying a vectorization transform to the matrix.

OK, so that's pretty straightforward, in fact you could have found all of that out from Wikipedia, so we're done right?

Actually, no we're not, in fact the above explanation, while correct, is a major cause of confusion.

For instance, you may be drawing the conclusion that a row major matrix is the transpose of a column matrix and you would be right, unfortunately you would also be wrong.

So why is that? Here are some of the reasons (in no particular order):

Row and column major only really describes array element indexing.
The elements of matrices and vectors are defined by their usage.
We have no definition of what differentiates a row and a column.
Vectors can be rows or columns.

To try and get a fuller understanding of row and column major, we need to take a quick detour through matrices and linear algebra:

Matrix multiplication is defined as a series of dot products between the rows of the matrix on the left and the columns of the matrix on the right with the resulting scalar values stored in the element with the row index from the left matrix and the column index from the right matrix. For instance the dot product of left row 2 and right column 3 ends up in the element at row 2, column 3.

Vector transformation uses the same rows on the left by columns on the right ordering as matrix multiplication.

While implementations may swap left and right parameters, from a purely mathematical point of view, row vectors must be the left side of equations and column vectors must be the right side of equations.

Now lets try a slightly different formulation of the row and column major definition which maps better to matrix maths:

A row major matrix is an array of row vectors and a column major matrix is an array of column vectors.

We also need to decide whether non-matrix vectors are going to be row vectors or column vectors. We could decide to be consistent and use row vectors with row major matrices and column vectors with column major matrices, but switching them can have benefits.

This leaves us with 4 potential scenarios, each with a slightly different set of ramifications:

Row vectors with row major matrices:

The vector is the left operand of the transform.
Vector transforms are dot products between the vector and the column elements of the matrix.
Matrix multiplication order is local space on the left and world space on the right.
Matrix vectors are the mapping of local space scalar units to world space vectors.

Row vectors with column major matrices:

The vector is the left operand of the transform.
Vector transforms are dot products between the vector and the vectors of the matrix.
Matrix multiplication order is local space on the left and world space on the right.
Matrix vectors are the planes of the local frame of reference.

Column vectors with row major matrices:

The vector is the right operand of the transform.
Vector transforms are dot products between the vectors of the matrix and the vector.
Matrix multiplication order is world space on the left and local space on the right.
Matrix vectors are the planes of the local frame of reference.

Column vectors with column major matrices:

The vector is the right operand of the transform.
Vector transforms are dot products between the row elements of the matrix and the vector.
Matrix multiplication order is world space on the left and local space on the right.
Matrix vectors are the mapping of local space scalar units to world space vectors.

The scenarios influence the memory layout of the matrices in the following combinations:

Scenarios 1 and 4 are identical.
Scenarios 2 and 3 are identical.
Scenarios 1 and 4 are transposed relative to scenarios 2 and 3.

So which scenario is best? Well, I'm on slightly dangerous ground here, but I'd say row vectors with column major matrices for the following reasons:

Vector matrix transforms can be performed with vector register aligned dot products.
Affine transforms (the most common type) can be represented as a matrix with only 3 column vectors, making the data compact while still aligned for vector registers.
A left to right, local to world transformation order matches the left to right reading order most familiar in the western world.
I view transformations in terms of planes, so having the matrix vectors be planes feels natural.

This is of course just my opinion and there are counter arguments which may favour other scenarios.

Ellipsoids

2011-04-10T14:35:00.000+01:00

I'm trying to collide 2 ellipsoids, anybody got any idea how to do this?

The ellipsoids are represented as 4x3 matrices and obviously the detection can be 'reduced' to an ellipsoid to unit sphere detection problem by transforming one ellipsoid by the inverse of the other, but that doesn't really help.

Colliding ellipsoids with points, lines, rays, planes, polygons is relatively trivial, but against spheres and ellipsoids seems fairly intractable to a linear algebra approach.

Ideally I'd like to include skewed ellipsoids.

A brief rant about splines

2011-03-27T16:14:00.002+01:00

Questions that should be asked if complaints are being made about spline processing performance:

If all your rational point weights are set to the same value, why are you using a rational spline?

If all your knot spacings are going to be the same, why are you using a non-uniform evaluator?

If your leading and trailing knot spacings are zero (to make the spline pass through the start and end control points) and all your other knot spacings are the same, why are you using a non-uniform evaluator? You could have just duplicated the start and end control points and saved both memory and processing.

If you need to evaluate the tangent to the spline, why are you actually evaluating points on the spline and taking the delta? This is slow and inaccurate.

If you need to evaluate the tangent to the spline and the position at the same t value, why are you not solving them with the same calculation? Evaluating them at the same time is little if any extra cost to just evaluating one of them.

Why are you using a recursive evaluation? Non-recursive is considerably faster regardless of the degree of the curve.

Why are you using an iterative evaluation if you are only ever evaluating cubic (or any other order known in advance) splines? Non-iterative is considerably faster than iterative, which itself is faster than recursive.

Why are you doing a linear search from the start of the knot vector to find your knot spans? A binary search is a good general purpose implementation and should be considerably quicker. If you happen to know that access will be sequential, a cached sequential search 'might' also be a good idea. If you are resorting to caching that last evaluation, something is wrong elsewhere in your code.

Are you dynamically changing your knot spacings? If not then consider pre-calculating span matrices.

For splines only (not surfaces): Are your control points and knot vector static? If so, then consider pre-calculating span matrices AND pre-multiplying these matrices by the control points, your spline evaluation cannot get much faster once you have done this.

Who's afraid of non-uniform rational b-splines (NURBS)?

2011-03-11T16:05:00.006Z

NURBS (non-uniform rational b-splines)

I remember being told by 3 separate artists that I must be incompetent for declining to do real-time NURBS on a pre-launch PS1 project (the PS1 didn't even have hardware floating point support). Despite this, I do have a fondness for NURBS.

These days, fast floating-point, hardware vector support and shaders are common and while NURBS may not be my first choice of spline or surface form, they are an option for real-time graphics. I would still probably avoid real-time collision directly against a NURBS surface and would definitely avoid NURBS based CSG in a real-time environment.

Before I start, I want to dispel a common misconception: while the definition of generalized NURBS may be recursive, the evaluation of a NURBS point does not have to be. The actual number of control points and knot values (see below) required to evaluate a point is only a fraction of the number that will generally be visited by a recursive evaluator.

First some definitions:

The degree of a curve is the highest power that the interpolation factor t is raised by.
The order of a curve is the degree plus one.
The minimum number of control points for a curve is the same as the order of the curve.
The number of elements in a knot vector is equal to the number of control points plus the degree of the curve minus one.
The evaluation of a point on a curve of order O requires O control points and 2*(O-1) knot values.
While NURBS can be of any degree and any dimensionality, they are typically degree 3 (cubic) with 3D control points.

Now let's take a look at NURBS and try to minimise the mathematical mumbo-jumbo. We'll start by breaking down the name NURBS into it's components:

NU (non-uniform):

In addition to an array of control points, the spline has an array of velocity control values known as a 'knot vector'. This knot vector consists of a series of increasing 'distance' values which control the speed at which each section of the spline is traversed. The gap between each pair of knot values is known as a 'span'. The knot vector can also be used to create discontinuities in the spline by having sections with zero distance between them.

Given a t value of 0 to 1, the section of the spline which needs to be evaluated can be found by multiplying t by the value of the last element of the knot vector and then scanning the list of knot vector values to find the index of the lowest distance value less than or equal to the scaled t value. A binary search is recommended to find this index value.

R (rational):

This means that each control point has a weight value associated with it.

When evaluating the spline, each control point is multiplied by it's weight value and then both the scaled control point and the weight are interpolated by the spline. Finally, the interpolated point is divided by the interpolated weight value.

The term 'rational' refers to the ratio between the weight scaled control point data and the weight. Rational control points can model shapes, such as a perfect circle, which cannot be produced with non-rational control points.

BS (b-spline):

Not to be confused with a Bezier curve, this is the real meat of the spline (pseudo-code for cubic NURBS can be found below).

Alright already, what about some pseudo-code for cubic NURBS evaluation?

Given:
  cp[]  : an array of control points (x, y, z components and w weight)
  kv[]  : an array of knot values
  bw[4] : storage for the blend weights produced by the evaluation
  n     : the number of control points
  t     : the spline interpolation value in the range 0 : kv[n+1]
  i     : the knot vector span index such that kv[i] <= t < kv[i+1]

Select active knot values:
 k0 = kv[ i-2 ]
 k1 = kv[ i-1 ]
 k2 = kv[ i ]
 k3 = kv[ i+1 ]
 k4 = kv[ i+2 ]
 k5 = kv[ i+3 ]

Select active control points:
  c0 = cp[ i-2 ]
  c1 = cp[ i-1 ]
  c2 = cp[ i ]
  c3 = cp[ i+1 ]

Calculate some common factors:
 l = ( ( k3 - t ) / ( ( k3 - k2 ) * ( k3 - k1 ) ) )
 r = ( ( t - k2 ) / ( ( k3 - k2 ) * ( k4 - k2 ) ) )
 a = ( ( r * ( t - k2 ) ) / ( k5 - k2 ) )
 b = ( ( ( r * ( k4 - t ) ) + ( l * ( t - k1 ) ) ) / ( k4 - k1 ) )
 c = ( ( l * ( k3 - t ) ) / ( k3 - k0 ) )

Calculate the blend weights:
 bw[ 0 ] = ( c * ( k3 - t ) )
 bw[ 1 ] = ( ( b * ( k4 - t ) ) + ( c * ( t - k0 ) ) )
 bw[ 2 ] = ( ( a * ( k5 - t ) ) + ( b * ( t - k1 ) ) )
 bw[ 3 ] = ( a * ( t - k2 ) )

Evaluate the blended control point weights:
  bw[ 0 ] *= c0.w
  bw[ 1 ] *= c1.w
  bw[ 2 ] *= c2.w
  bw[ 3 ] *= c3.w
  w = bw[ 0 ] + bw[ 1 ] + bw[ 2 ] + bw[ 3 ]

Evaluate the final point:
  x = (( c0.x*bw[0] )+( c1.x*bw[1] )+( c2.x*bw[2] )+( c3.x*bw[3] ))/w
  y = (( c0.y*bw[0] )+( c1.y*bw[1] )+( c2.y*bw[2] )+( c3.y*bw[3] ))/w
  z = (( c0.z*bw[0] )+( c1.z*bw[1] )+( c2.z*bw[2] )+( c3.z*bw[3] ))/w

The matrix form of this can be found by extracting 't' from the formula, giving 4 vector4 constants (one for each power of t from 0 to 3).

By the way, evaluating the partial first derivative (and tangent) can share a considerable amount of the work required to evaluate the point, if a reader expresses an interest in this I will add that formulation.

NURBS lend themselves to a considerable amount of pre-calculation:

Assuming that the knot vector is not being modified, each knot can be pre-evaluated to a matrix.

Even if the knot vector is being modified (which is unlikely in the context of game rendering), recalculating the knot matrices is low cost, though if large sections of the knot vector are frequently being modified it is probably preferable to stick to a direct rather than pre-calculated spline evaluation.

In this scenario, evaluating a NURBS point will require 2 scalar multiplies, 1 scalar reciprocal, 1 vector3 scale and 5 vector4 by matrix44 multiplies.

Modifying the control points does not require recalculation of the knot evaluation matrices.

Assuming that both the control points AND the knot vector are not being modified, the control points can be potentially be pre-multiplied by the knot evaluation matrices.

For this to be possible, the degree of the curve must be the same as the dimensions of the control points, for example a cubic (degree 3) NURBS with 3D control points or a quadratic (degree 2) NURBS with 2D control points.

In this ideal scenario, evaluating a NURBS point will require 2 scalar multiplies, 1 scalar reciprocal, 1 vector3 scale and 1 vector4 by matrix44 multiply.

While it might appear that pre-multiplying the control points by their weights is an optimization, in general it is not.

In the dynamic calculation shown above, it can be seen that the multiply is absorbed into the general weight evaluation. In a pre-calculated scenario, the weights are best multiplied into the knot evaluation matrix.

Considerations for Stereoscopic 3D

2011-02-15T15:51:00.001Z

There are a lot of factors which influence successful 3D stereoscopy, here are the ones I consider most important:

Physical dimensions of the screen
Physical position of the observer
Physical eye separation (interocular distance)
Physical eye rotation, required field of view, zoom and scene depth
Physical eye focus and speed of depth transitions
Scene clipping and depth
Lighting

Physical screen dimensions and the observer position:

Together the dimensions of the screen and position of the observer give the neutral (unforced) field of view.

Neutral horizontal FOV = 2 * atan( screen width / ( observer distance from the screen * 2 ) )

Typically, for a comfortable viewing position, the FOV is about 55 degrees (or less). This means that the width of the view-port and the distance of the eye point from the view plane should be roughly the same for a neutral projection.

For the screen to be perceived as a window, the ratios of screen width, screen height and observer position need to be recreated in the projection.

Eye separation:

The human interocular distance (eye separation) is generally assumed to be about 2.5inches or 6.5cm and this should be factored into the projection setup using the same ratio that was used to map the screen to the view-port. This should result in an eye separation value which is approximately 1/6th of the width of the view-port.

Eye rotation, required field of view, zoom and scene depth:

The limits of human eye rotation (both convergence and divergence) have a large impact on the ability to resolve a 3D view. Most people can converge their eyes (cross them) to a significant degree, however very few can diverge them beyond the parallel. The implications of this are:

Moving the view plane closer to the eye (or increasing the FOV or view-port size):

Increases the perception of depth for objects behind the view plane but makes the scene harder to resolve.
Increases the apparent speed of object depth changes which can be hard to resolve.
Exponentially increases clipping issues for objects infront of the screen.
Increases the eye divergence required to resolve objects behind the view plane.
Reduces the resolvable scene depth (too much depth and the 3D illusion will fail).
Reduces the amount of depth available and resolvable infront of the view plane.

Moving the view plane further from the eye (or decreasing the FOV of view-port size):

Increases the depth available to make objects 'hang' infront of the screen.
Reduces clipping issues for objects infront of the screen.
Reduces the apparent speed of object depth changes and so is easier to resolve.
Makes the scene generally easier to resolve, but reduces depth perception for objects behind the view plane.

In general, view planes closer than the neutral depth should only be used in enclosed areas with limited depth to avoid breaking the stereoscopic illusion.

Eye focus and speed of depth transitions:

Eye focus has only an indirect influence on resolving stereoscopic images, but it is still an important consideration.

Stereoscopic 3D requires that the eyes focus remain on the surface of the screen independant of the eye rotation. HOWEVER, normally the eye focus and rotation are coordinated and separating the two requires some conscious effort. As the disparity between the eye focus distance and rotation distance increased, so does the likelyhood of eye-strain, headaches and failure to resolve the image.

Fast object depth transitions should be avoided as the eyes need time to adjust their rotation while maintaining their focus on the screen.

The amount of time that important scene objects are not on or near the view plane should be limited so that the eyes can have time to 'relax'.

Scene clipping:

Objects which only appear to one eye can break the 3D illusion. This is not a significant issue for objects behind the view plane, but can be a serious problem for objects infront of the screen.

I recommend that a a common object rejection clipping frustum be used for both the left and right scenes. The rejection frustum should have its base directly between the eyes and include the maximum far extents of both eye clip frustums. This is scene clipping frustum not to be confused with the polygon clipping frustums produced by the stereoscopic projection matrices.

Lighting:

The position of the left and right eyes should be used for rendering the lighting of the left and right scenes respectively. Using only a single eye position will work reasonably for objects behind the view plane, it will break (or at least strain) the 3D illusion for objects infront of the view plane.

In summary:

To reduce the possibility of eye-strain and improve the experience of 3D:

The speed at which objects transition between depths should be controlled to let the eyes adjust so that the object can be continuously resolved.
The amount of time that objects are infront of the view plane should be inversely proportional to their distance infront of the view plane. This reduces eye-strain and helps preserve the integrity of the resolution of the overall scene.
Don't move objects too close to the observer.
A narrower field of view is usually (but not always) preferable to a wide field of view.
Most action should take place at or near the view plane, ideally just behind it.

Ultimately, the golden rule for successful stereoscopic 3D is: use it sparingly and make it count. Constantly having objects hanging infront of the screen reduces the WOW factor compared to the occasional targetted use of objects coming out of the screen.

Derivation of 'the ultimate projection calculation'

2011-01-30T22:03:00.004Z

For the morbidly curious (among whom I count myself) here is the deriviation of 'the ultimate projection calculation'.

By the way, a prettier version (with diagrams!) of this derivation can be found on Martin Ridgers's blog here.

First a reminder of our variables:

  s   : stereo-scopic 3D eye separation
  p   : perspective (0 == parallel, 1 == perspective)
  n, f: near and far z clip depths
  w, h: width and height of viewport (at depth v_z)
  v_xyz   : center of viewport
  e_xyz   : center of projection (eye position)

Let's start with the X and Y planes:

We know that the X and Y planes are basically going to be the same except that the parameters will be x and y and the incoming x and y values will map to the x and y components of the X and Y planes. Given that, we only really need to derive the X or the Y plane because the other plane will just be a remapping. So, let's derive Y:

First we'll apply a skew so that the output y values for the points v and e are both zero. Given that we are talking about a zero value, we don't need to consider the scaling due to the perspective divide by w, after all, 0/w is still zero. The Y transformation to do this will be of the form:

  v_xY_x + v_yY_y + v_zY_z + Y_w = 0
  e_xY_x + e_yY_y + e_zY_z + Y_w = 0

This has too many variables to solve, however at this point given that we know that x has no influence on y (at least in this formulation), Y_x should be 0. We can also assume that Y_y is 1 so that we are simply passing the input y value through the equation without any scaling. This leaves us with:

  v_y + v_zY_z + Y_w = 0
  e_y + e_zY_z + Y_w = 0

This is clearly a simultaneous equation (two equations with two unknowns) and can be simply solved to:

  Y_z = ( e_y - v_y ) / ( v_z - e_z )
  Y_w = ( v_ye_z - e_yv_z ) / ( v_z - e_z )

So, now we have derived the skew contribution, we need to scale the viewport into the -1:+1 range. We want v_y + h/2 to map to +1 (the top edge of our frustum) and v_y - h/2 to map to -1 (the bottom edge of our frustum). As we have not applied any scaling to our Y plane at this point and we know that v_y maps to 0, the top edge of our viewport currently maps to +h/2 and the bottom edge to -h/2. Simply dividing our existing Y plane by h/2 will correctly map them to the -1:+1 range. This modifies our Y plane formula to:

  Y_y = 2 / h
  Y_z = 2( e_y - v_y ) / h( v_z - e_z )
  Y_w = 2( v_ye_z - e_yv_z ) / h( v_z - e_z )

This leaves us with the problem of the perspective divide. Our nice -1:+1 range is going to be divided by the w value resulting from the W plane. Given that we can decide what the w value will be at the viewport depth, let's just remember to make it 1 when we come to our W plane formula and leave our Y and X plane formulas unchanged.

By the way, it should be clear from this that simply adding and subtracting half the stereo 3D eye seperation 's' from e_x value in the X plane formula will give us left and right eye versions of our X plane as follows:

  X_x = 2 / w
  X_z = ( 2( e_x - v_x ) + ( ±s ) ) / w( v_z - e_z )
  X_w = ( 2( v_xe_z - e_xv_z ) - ( ±s )v_z ) / w( v_z - e_z )

Let's move on to the W plane:

The W plane controls perspective, so we are going to need to take account of our perspective control parameter 'p'. We are also going to have to ensure that the w value on the viewport plane is 1 as required by our formula for the X and Y planes.

A parallel projection W plane is going to look like:

  W_x = 0
  W_y = 0
  W_z = 0
  W_w = 1

This will have no perspective and the w value at the viewport depth will be 1 (let's face it, the w value will be 1 at any depth.

A perspective W plane needs to pass through the eye point at depth e_z so we can start with a normalized plane of the form:

  W_x = 0
  W_y = 0
  W_z = 1
  W_w = -e_z

This will offset z relative to the eye depth, but the w value at the viewport depth will be v_z * 1 - e_z, not the 1 that we require. This can be remedied by simply dividing the perspective W plane by ( v_z - e_z ) as follows:

  W_x= 0
  W_y= 0
  W_z= 1 / ( v_z - e_z )
  W_w= -e_z / ( v_z - e_z )

Now we have a perspective projection W plane with a w value of 1 at the viewport depth and a parallel projection W plane also with a w value of 1 at the viewport depth. Our perspective control parameter 'p' can simply interpolate between these 2 planes, giving us a W plane formula of:

  W_x = 0
  W_y = 0
  W_z = p / ( v_z - e_z )
  W_w = ( v_z(1-p) - e_z ) / ( v_z - e_z )

Now we're only left with the Z plane:

We know that x and y will have no contribution to our Z plane and so we are only left to solve Z_z and Z_w. This is going to be slightly more complicated than the solutions for the other planes as we need to map near (n) and far (f) z values to -1 and +1 respectively after perspective divides by different w values, rather than by 1 as we were able to use with the X and Y planes. This means that we have to solve:

  ( fZ_z + Z_w ) / ( fW_z + W_w ) = +1
  ( nZ_z + Z_w ) / ( nW_z + W_w ) = -1

Thankfully this is just another simultaneous equation and can be trivially solved to:

  Z_z = ( 2W_w + W_z( f + n ) ) / ( f - n )
  Z_w = -( W_w( f + n ) + 2fnW_z ) / ( f - n )

Which then expands (by substituting the formulas for W_z and W_w and rearranging) to a Z plane equation of:

  Z_x = 0
  Z_y = 0
  Z_z = ( 2( v_z(1-p) - e_z ) + p( f + n ) ) / ( ( f - n )( v_z - e_z ) )
  Z_w = -( ( v_z(1-p) - e_z )( f + n ) + 2fnp ) / ( ( f - n )( v_z - e_z ) )

Anything else?

Not a lot really except that if you want the w value at the viewport depth to be other than 1, multiply the entire matrix by the new value.

If you want to make W_z == 1 (as required by some SDKs and hardware), divide the entire matrix by W_z, though this will not work for parallel projection matrices and will have the side effect of modifying the w value at the viewport depth.

Frustum plane extraction

2011-01-29T15:08:00.000Z

The following extracts the frustum face planes from a unit (ogl) matrix.

Left face plane:

    L_x = W_x + X_x
    L_y = W_y + X_y
    L_z = W_z + X_z
    L_w = W_w + X_w

Right face plane:

    R_x = W_x - X_x
    R_y = W_y - X_y
    R_z = W_z - X_z
    R_w = W_w - X_w

Bottom face plane:

    B_x = W_x + Y_x
    B_y = W_y + Y_y
    B_z = W_z + Y_z
    B_w = W_w + Y_w

Top face plane:

    T_x = W_x - Y_x
    T_y = W_y - Y_y
    T_z = W_z - Y_z
    T_w = W_w - Y_w

Near face plane:

    N_x = W_x + Z_x
    N_y = W_y + Z_y
    N_z = W_z + Z_z
    N_w = W_w + Z_w

Far face plane:

    F_x = W_x - Z_x
    F_y = W_y - Z_y
    F_z = W_z - Z_z
    F_w = W_w - Z_w

This also works for DX matrices except for the near plane which should be as shown below.

Near face plane (DX):

    N_x = Z_x
    N_y = Z_y
    N_z = Z_z
    N_w = Z_w

Notes:

Obviously this formula can be used to extract the frustum faces from any sub-frustum by first applying a frustum remap to the matrix.

Extracting the vertices of the frustum is usually best done by inverting the frustum matrix and summing the planes to simulate transforming a unit (-1:+1) cube.

The frustum vertices can also be extracted by intersection of the frustum planes and a optimized vertex extractor of this form can have comparable performance to the inversion and summing method.

My favourite spline

2011-01-24T14:14:00.004Z

For a quick diversion, I thought I'd mention my favourite spline which is a variant on a non-uniform Catmull-Rom spline.

A Catmull-Rom spline is commonly implemented as a uniform cubic Hermite curve with tangents:

    U=(P_i+1-P_i-1)/6 and V=(P_i+2-P_i)/6

Given a knot vector K with knot spacing equal to the distance between the control points, replacing the tangents with:

    U=(P_i+1-P_i-1)(K_i+1-K_i)/3(K_i+1-K_i-1)
    V=(P_i+2-P_i)(K_i+1-K_i)/3(K_i+2-K_i)

Produces a very nice non-uniform interpolating spline with very little acceleration because the tangents take account of the accelleration due to the interpolation of the knot vector.

Frustum remapping

2011-01-22T17:09:00.003Z

So, we have a unit frustum but we want to remap it because:

We want to remap it to a non-unit frustum (like DirectX).
We want to use a sub-section of it for tiled rendering.
We want to invert depth for better floating point depth buffer accuracy.
We want the rendering to be upside down.
...

Remapping the frustum is simply a matter of combining various ratios of the W plane and the X, Y and Z planes to shift the center point and rescale the frustum.

Given 2 aabb's s and t representing the required remapping from source to target, the formula is as follows:

Calculate the ratios:

    sx = 1.0 / ( s.max.x - s.min.x )
    sy = 1.0 / ( s.max.y - s.min.y )
    sz = 1.0 / ( s.max.z - s.min.z )
    xx = ( t.max.x - t.min.x ) * sx
    yy = ( t.max.y - t.min.y ) * sy
    zz = ( t.max.z - t.min.z ) * sz
    xw = ( s.max.x * t.min.x - t.max.x * s.min.x ) * sx
    yw = ( s.max.y * t.min.y - t.max.y * s.min.y ) * sy
    zw = ( s.max.z * t.min.z - t.max.z * s.min.z ) * sz

Apply the remapping:

    X_x = X_x * xx + W_x * xw
    X_y = X_y * xx + W_y * xw
    X_z = X_z * xx + W_z * xw
    X_w = X_w * xx + W_w * xw
    Y_x = Y_x * yy + W_x * yw
    Y_y = Y_y * yy + W_y * yw
    Y_z = Y_z * yy + W_z * yw
    Y_w = Y_w * yy + W_w * yw
    Z_x = Z_x * zz + W_x * zw
    Z_y = Z_y * zz + W_y * zw
    Z_z = Z_z * zz + W_z * zw
    Z_w = Z_w * zz + W_w * zw

How about some examples?

To map a unit (ogl) frustum to a DX frustum:

      s.min = { -1, -1, -1 }
      s.max = {  1,  1,  1 }
      t.min = { -1, -1,  0 }
      t.max = {  1,  1,  1 }

To map a DX frustum to a unit (ogl) frustum:

      s.min = { -1, -1,  0 }
      s.max = {  1,  1,  1 }
      t.min = { -1, -1, -1 }
      t.max = {  1,  1,  1 }

To mirror a unit frustum (inverts polygon winding order):

      s.min = {  1, -1, -1 }
      s.max = { -1,  1,  1 }
      t.min = { -1, -1, -1 }
      t.max = {  1,  1,  1 }

To make a unit frustum upside down (inverts polygon winding order):

      s.min = { -1,  1, -1 }
      s.max = {  1, -1,  1 }
      t.min = { -1, -1, -1 }
      t.max = {  1,  1,  1 }

To invert the depth (but not the perspective or handedness) a unit frustum:

      s.min = { -1, -1,  1 }
      s.max = {  1,  1, -1 }
      t.min = { -1, -1, -1 }
      t.max = {  1,  1,  1 }

To map the far, top, right octant of a unit frustum to a DX frustum:

      s.min = {  0,  0,  0 }
      s.max = {  1,  1,  1 }
      t.min = { -1, -1,  0 }
      t.max = {  1,  1,  1 }

To map the near, bottom, left octant of a DX frustum to a unit frustum:

      s.min = { -1, -1,  0   }
      s.max = {  0,  0,  0.5 }
      t.min = { -1, -1, -1   }
      t.max = {  1,  1,  1   }

Fair enough, but where do those ratios come from?

The easiest way of understanding this is with a unit frustum. There are 2 parts to the remapping, translating the frustum origin and scaling the frustum. I'll only cover the x axis, but the same logic applies to the y and z axes.

Translation:

After the perspective divide by w, the left edge of a unit frustum maps to x = -1 and the right to x = +1. This means that before the perspective divide, the left edge x value was -w and the right was +w. From this, you should be able to see that the W plane basically defines a unit of distance from the center of a unit frustum.

So, to move the center line of the frustum to the left edge of the frustum, simply subtract the W plane from the X plane, everything which previously had a x value of 0 after the perspective devide will now have a value of -1.

That should be enough to get you started with moving the origin of the frustum.

Scaling:

Scaling the frustum is even simpler. To perform a 2x horizontal zoom, we simply need to double the X plane, this will double all the output x values, mapping values of 0.5 to 1.0.

So why change my projection functions?

2011-01-22T15:31:00.003Z

This article refers to the projection matrix formulation in this post.

Why is this any better than what I currently use?

It is easy to visualize and manipulate.
It has a lot more functionality.
It can be interpolated.
It supports stereoscopic 3D.
It simplifies shadow mapping adjustment.
It correctly supports the skews required for oblique projections.
It simplifies the use of offset projections (which are generally under-used because they are not understood).
It can reduce the number of code paths because the same calculation can be used for multiple matrix forms.
...

OK, any specifics?

Camera pull:
Dolly zoom (aka vertigo zoom or tracking zoom):
Planar zooms:
Simple zoom multiples:
Blending between parallel and perspective projections:
Oblique projections:
Controlling stereoscopic focus:
Driving user attention:

But what about FOV and the other controls I'm used to?

These can all be implemented as helper functions to set the parameters this form requires. A modifier queue should already form part of most good camera systems and all your existing controls should be easy to adapt to this form.

Why only 'almost'?

2011-01-22T13:18:00.002Z

OK, so why is the formulation given in my previous post only 'almost' the ultimate projection matrix? Well, it does not account for the following:

Non-unit frustums (i.e. x, y, and z extents other than -1:+1).
Handedness (the projection is left handed).
Depth biasing.
Infinite projections.
Skewed near and far clipping planes.

But why not?

Non-unit frustums (i.e. x, y, and z extents other than -1:+1):

A unit frustum is the most neutral formulation.
Even if a unit frustum is not the final desired projection, it is a useful intermediate matrix for use with collision detection.
The remapping from a unit frustum to another set of bounds (such as the DirectX form of x[-1:+1, y[-1:+1], z[0:+1]) is a trivial modification which can be applied after this formulation. Having a seperate remap function is useful as part of a seperate pre-viewport transform object and can be used to create sub-frustums or invert bounds directions. An explanation of how the remapping is performed will be given in an upcoming post.

Handedness (the projection is left handed):

Does not generally contribute to creative camera control.
Cannot be interpolated.
To convert the matrix to a right handed projection, simply negate the z components of each plane of the matrix.

Depth biasing:

Is a special case which only modifies the Z plane.
Does not generally contribute to creative camera control.
There are alternative methods of applying depth biasing.
Could be added to the formulation, but:

Is a geometry and rendering problem rather than a camera control.
Requires a degree of knowledge of the available depth precision.

Infinite projections:

Are a special case which only modifies the Z plane.
Do not generally contribute to creative camera control.
Can be useful for shadow volumes and similar, but:

Require care with precision and numerical stability.
Cause issues with functions expecting a finite projection.
Are best handled as special cases or matrix modifiers.

Skewed near and far clipping planes:

Is a special case which only modifies the Z plane.
Do not generally contribute to creative camera control.
Can be useful for portal rendering and similar, but:

Require decisions about which plane to modify and how to handle the remaining depth clip plane.
Are hard to visualize due to the relationship between the near and far plane.
Do not produce a depth buffer consistent with the rest of the scene.
There are alternative methods of dealing with portals.
Are best handled as special cases or matrix modifiers.

If you want to know more about matrix biasing, infinite projections and skewed near and far clipping planes, I suggest you check out the excellent explanation by Eric Lengyel which can be found here.

Up next: OK fine, so why is this any better than what I currently use?

The ultimate projection calculation (almost)

2010-11-25T23:27:00.008Z

The following builds perspective and parallel projection matrices with support for center offsets and stereoscopic 3D.

Variables:

    s   : stereo-scopic 3D eye separation
    p   : perspective (0 == parallel, 1 == perspective)
    n, f: near and far z clip depths
    w, h: width and height of viewport (at depth v_z)
    v_xyz   : center of viewport
    e_xyz   : center of projection (eye position)

W-Plane:

    W_x = 0
    W_y = 0
    W_z = p/(v_z-e_z)
    W_w = (v_z(1-p)-e_z)/(v_z-e_z)

Z-Plane:

    Z_x = 0
    Z_y = 0
    Z_z = (2(v_z(1-p)-e_z)+p(f+n))/((f-n)(v_z-e_z))
    Z_w = -((v_z(1-p)-e_z)(f+n)+2fnp)/((f-n)(v_z-e_z))

Y-Plane:

    Y_x = 0
    Y_y = 2/h
    Y_z = 2(e_y-v_y)/(h(v_z-e_z))
    Y_w = 2(v_ye_z-e_yv_z)/(h(v_z-e_z))

X-Plane:

    X_x = 2/w
    X_y = 0
    X_z = (2(e_x-v_x)+[±s])/(w(v_z-e_z))
    X_w = (2(v_xe_z-e_xv_z)-[±sv_z])/(w(v_z-e_z))

Notes:

The projection produced by this formula has x, y and z extents of -1:+1.

The perspective control value p is not restricted to integer values.

The view plane is defined by z. Objects on the view plane will have a homogeneous w value of 1.0 after the transform.

Upcoming posts will explain what this formulation does, what it does not do and why you should use it.

What it's all about

2010-11-13T20:18:00.000Z

I'm a games programmer. I've been a games programmer for a very long long time. As of a day or two ago, I'm "between jobs".

So what is this blog for? Well, it's a place to gather my thoughts and hopefully share some of what I have learned and what I find interesting.

My first subject is going to be one that has frustrated me for some time - virtual camera projections. Most graphics programmers will have used projection matrices, but few know much about the full range of features that can be captured by this matrix, not just to simplify common rendering problems but to grab and drive the attention of the viewer. Common implementations of the projection matrix supplied by SDKs and off the shelf implementations often attempt to simplify the matrix construction, but that very simplification hides the power of this matrix and frequently results in more complicated code and multiple matrix construction functions.

The aim of my first discussion is to show the possibilities embedded within the projection matrix.

First post

2010-10-14T03:05:00.000+01:00

Details to follow...

And what happened?

Fast 2D and 3D Hilbert curves and Morton codes

Row major and column major matrices

Row major and column major matrices (or proving that black is white)

Ellipsoids

A brief rant about splines

Who's afraid of non-uniform rational b-splines (NURBS)?

NURBS (non-uniform rational b-splines)

NURBS lend themselves to a considerable amount of pre-calculation:

Considerations for Stereoscopic 3D

Physical screen dimensions and the observer position:

Eye separation:

Eye rotation, required field of view, zoom and scene depth:

Eye focus and speed of depth transitions:

Scene clipping:

Lighting:

In summary:

Derivation of 'the ultimate projection calculation'

First a reminder of our variables:

Let's start with the X and Y planes:

Let's move on to the W plane:

Now we're only left with the Z plane:

Anything else?

Frustum plane extraction

Left face plane:

Right face plane:

Bottom face plane:

Top face plane:

Near face plane:

Far face plane:

Near face plane (DX):

Notes:

My favourite spline

Frustum remapping

So why change my projection functions?

This article refers to the projection matrix formulation in this post.

Why is this any better than what I currently use?

OK, any specifics?

Camera pull:

Dolly zoom (aka vertigo zoom or tracking zoom):

Planar zooms:

Simple zoom multiples:

Blending between parallel and perspective projections:

Oblique projections:

Controlling stereoscopic focus:

Driving user attention:

But what about FOV and the other controls I'm used to?

Why only 'almost'?

The ultimate projection calculation (almost)

Variables:

W-Plane:

Z-Plane:

Y-Plane:

X-Plane:

Notes:

What it's all about

First post