OpenGL matrix setup for tiled rendering

Question

After reading datenwolf's 2011 answer concerning tile-based render setup in OpenGL, I attempted to implement his solution. The source image looks like this (at 800 x 600)

The resulting image with 2x2 tiles, each tile at 800 x 600 per tile looks like this.

As you can see they don't exactly match, though I can see something vaguely interesting has happened. I'm sure I've made an elementary error somewhere but I can't quite see it.

I'm doing 4 passes where:

w, h are 2,2 (2x2 tiles)    
x, y are (0,0) (1,0) (0,1) and (1,1) in each of the 4 passes    
MyFov is 1.30899692 (75 degrees)
MyWindowWidth, MyWindowHeight are 800, 600
MyNearPlane, MyFarPlane are 0.1, 200.0

The algorithm to calculate the frustum for each tile is:

auto aspect = static_cast<float>(MyWindowWidth) / static_cast<float>(MyWindowHeight);
auto right = -0.5f * Math::Tan(MyFov) * MyShaderData.Camera_NearPlane;
auto left = -right;
auto top = aspect * right;
auto bottom = -top;
auto shift_X = (right - left) / static_cast<float>(w);
auto shift_Y = (top - bottom) / static_cast<float>(h);
auto frustum = Math::Frustum(left   + shift_X * static_cast<float>(x),
                             left   + shift_X * static_cast<float>(x + 1), 
                             bottom + shift_Y * static_cast<float>(y),
                             bottom + shift_Y * static_cast<float>(y + 1),
                             MyShaderData.Camera_NearPlane,
                             MyShaderData.Camera_FarPlane);

, where Math::Frustum is:

template<class T>
Matrix4x4<T> Frustum(T left, T right, T bottom, T top, T nearPlane, T farPlane)
{
    Matrix4x4<T> r(InitialiseAs::InitialiseZero);

    r.m11 = (static_cast<T>(2) * nearPlane) / (right - left);
    r.m22 = (static_cast<T>(2) * nearPlane) / (top - bottom);
    r.m31 = (right + left) / (right - left);
    r.m32 = (top + bottom) / (top - bottom);
    r.m33 = -(farPlane + nearPlane) / (farPlane - nearPlane);
    r.m34 =   static_cast<T>(-1);
    r.m43 = -(static_cast<T>(2) * farPlane * nearPlane) / (farPlane - nearPlane);

    return r;
}

For completeness, my Matrx4x4 layout is:

struct
{
    T m11, m12, m13, m14;
    T m21, m22, m23, m24;
    T m31, m32, m33, m34;
    T m41, m42, m43, m44;
};

Can anyone spot my error?

Edit:

So derhass explained it to me - a much easier way of doing things is to simply scale and translate the projection matrix. For testing I modified my translation matrix, scaled up by 2x, as follows (changing translate for each tile):

auto scale = Math::Scale(2.f, 2.f, 1.f);
auto translate = Math::Translate(0.5f, 0.5f, 0.f);
auto projection = Math::Perspective(MyFov,
                                    static_cast<float>(MyWindowWidth) / static_cast<float>(MyWindowHeight),
                                    MyShaderData.Camera_NearPlane,
                                    MyShaderData.Camera_FarPlane);          

MyShaderData.Camera_Projection = scale * translate * projection;

The resulting image is below (stitching 4 together) - the discontinuities in the image are caused by the post processing I think, so that's another issue I might have to deal with at some point.

Answer 1

This isn't a real answer for the question, but it might be an useful alternative approach to what you are trying to solve here. In my opinion, datenwolf's solution in his answer to the stackoverflow question you are referring to is more comlicated that it needs to be. So I'm presenting my alternative here.

Forword: I assume standard OpenGL matrix conventions, so that the vertex transformation with matrix M is done as v'= M *v (like the fixed-function pipeline did).

When a scene is rendered with some projection matrix P, you can extract any axis-aligned sub-rectangle of said scene by applying a scale and transformation operation after the projection matrix is applied.

The key point is that the viewing volume is defined as the [-1,1]^3 cube in NDC space. The clip space (which is what P transforms the data to) is just the homogenous represenation of that volume. As the typical 4x4 transformation matrices are all working in homogenous space, we don't really need to care about w at all and simply can define the transformations as if we were in NDC space.

Since you only need some 2D tiling, z should be left as-is, and only some scale and translation in x and y is required. When composing transformations A and B into a single Matrix C as C=A*B, following the aforementioned conventions this results in B being applied first, and A last (since C*v == A*B*v == A*(B*v)). So to modify the result after projection, we have to pre-multiply some transformations to P and we are done:

P'=S(sx,sy,1) * T(tx,ty,0) * P

The construction of P' will work with any valid projection matrix P, no matter if it is a perspective or ortho transform. In the ortho case, what this does is quite clear. In the perspective case, this actually modifies both the field of view and also shifts the frustum to an asymmetric one.

When you want to tile the image into a grid of m times n segments. it is clear that sx=m and sy=n. As I did use the S * T order (by choice), T is applied before the scale, so for each tile, (tx,ty) is just the vector moving the center of the tile to the new center (which will be the origin). As NDC space is 2 units wide and tall, for a tile x,y, the transformation is

tx= - (-1 + 2/(2*m) + (2/m) * x)
ty= - (-1 + 2/(2*n) + (2/n) * y)
//     ^     ^         ^    
//     |     |         |
//     |     |         +- size of of each tile in NDC space
//     |     |
//     |     +- half the size (as the center offset)
//     |
//     +- left/bottom border of NDC space