Shading a concave cube as a convex cube with forced perspective

Question

Context

I'm building a graphics pipeline for voxel volumes. I'm using an existing game engine (Bevy) which provides a way to put an object in 3D space. In my application, the voxel volumes can be oriented arbitrarily and are not necessarily axis aligned with respect to the world space. The pipeline does front-face culling of these voxel volumes to allow the camera to go inside the bounds of the volume without clipping (imagine a player walking into a sparsely populated volume and looking around). The voxel volumes are rendered using ray marching, starting from a ray's contact point with the front face of the voxel volume.

The Problem

Now of course, this works if there are front faces, but as I said I'm culling them. So what really happens is my frag shader is shading the near side of the cube's back faces. What I need to do is: given the point of contact with the cube's back face, and the direction to the camera, find the point that would have been hit on the front face. This will allow the rest of the program to ray march as though the ray hit a "convex" cube. I've spent quite a while looking at this answer but so far I've been unable to implement it into my shader -- I think it's relevant: https://stackoverflow.com/questions/4248090/finding-the-length-of-a-ray-within-a-cube

Assumptions:

• Voxel volume axes are perpendicular to each other, but are not necessarily the same length
• A ray from the camera to a back face can hit any one of the 3 back faces, and will only pass through any one of the 3 front faces

Expected Result

Here's an example of what I expect. Below (first image) is what is rendered without "simulating" the front face positions. The second image is what would be rendered with correct front face simulation.

Partial Code

I've removed the irrelevant code. Right now what is rendered for a voxel_volume_size = (16, 16, 16) is 3 16x16 walls along the axes. Of course what should be seen is a 16x16x16 voxel cube.

voxel.vs

#version 450

layout(location = 0) in vec3 Vertex_Position;
layout(location = 1) in vec3 Vertex_Normal;
layout(location = 2) in vec2 Vertex_Uv;

layout(location = 0) out vec3 v_Position;
layout(location = 1) out vec3 v_Normal;
layout(location = 2) out vec2 v_Uv;

layout(set = 0, binding = 0) uniform Camera {
    mat4 ViewProj;
    mat4 View;
};

layout(set = 1, binding = 0) uniform Transform {
    mat4 Model;
};

void main() {
    v_Normal = Vertex_Normal;
    v_Position = Vertex_Position;
    v_Uv = Vertex_Uv;
    gl_Position = ViewProj * vec4((Model * vec4(Vertex_Position, 1.0)).xyz, 1.0);
}

voxel.fs

#version 450

layout(location = 0) in vec3 v_Position;
layout(location = 1) in vec3 v_Normal;
layout(location = 2) in vec3 v_Uv;

layout(location = 0) out vec4 o_Target;

layout(set = 0, binding = 0) uniform Camera {
    mat4 ViewProj;
    mat4 View;
};

layout(set = 1, binding = 0) uniform Transform {
    mat4 Model;
};

layout(set = 3, binding = 0) buffer VoxelVolume {
    vec3 voxel_volume_size;
};

void main(void) {
    vec3 Normal = mat3(Model) * v_Normal;

    vec3 scale = voxel_volume_size / 16.0;

    mat4 InverseView = inverse(View);
    vec3 CameraPosition = (Model * vec4(vec3(InverseView[3]), 0.)).xyz;

    vec3 BackFacePosition = v_Position;
    vec3 BackFaceModelPosition = (Model * vec4(BackFacePosition, 1.0)).xyz;
    vec3 BackFaceRayDirection = normalize(BackFaceModelPosition - CameraPosition);

    vec3 FrontFacePosition = ?;
    vec3 FrontFaceRayDirection = ?; // Differs from BackFaceRayDirection when front face is perpendicular to back face

    vec3 ScaledPosition = ((FrontFacePosition + (scale / 2.0)) / scale) * voxel_volume_size;

    o_Target = vec4(floor(ScaledPosition) / voxel_volume_size, 1.0);
}

I would be extremely grateful if someone could help me fill in those missing variables! Thanks.

Answer 1

I ended up solving this with a little of guesswork. I suppose this can be optimized further.

---

Matrix setup

First, some mental health. I declared a convert function:

vec3 convert(mat4 matrix, vec3 vector)
{
    return (matrix * vec4(vector, 1.0)).xyz;
}

And then in the fragment shader I have my matrices like this:

    mat4 ModelToWorld = Model;
    mat4 WorldToModel = inverse(Model);
    mat4 CameraToWorld = inverse(View);
    mat4 WorldToCamera = View;
    mat4 ModelToCamera = WorldToCamera * ModelToWorld;
    mat4 CameraToModel = WorldToModel * CameraToWorld;

I will prefix the variables with the space in which they are. For example:

    vec3 Model_BackFacePosition = v_Position;

And if I want that in Camera space I do this:

    vec3 Camera_BackFacePosition = convert(ModelToCamera, Model_BackFacePosition);

This is somewhat verbose, but it is less error prone. I can see if I'm converting from the wrong space using the wrong matrix. Or if I'm mixing vectors of different spaces unintentionally.

Makes sense? Ok.

---

Solution

We know:

    vec3 Model_BackFacePosition = v_Position;

That the is position, in model coordinates, we got from the vertex shader.

We also know that:

    vec3 Camera_RayOrigin = vec3(0.0, 0.0, 0.0);

Because the origin of the ray is the camera. And the position of the camera respect to the camera is 0. Let us pass that to Model space:

    vec3 Model_RayOrigin = convert(CameraToModel, Camera_RayOrigin);

And we also need the direction of the ray, which we get like this:

   vec3 Model_RayDirection = normalize(Model_BackFacePosition - Model_RayOrigin);

---

Now we need to figure out the normals of the front faces. There are at most three front faces (at least one). And they are aligned to the axis in Model space. So let us get them in Model space.

    vec3 Model_XN = vec3(-sign(Model_RayOrigin.x), 0.0, 0.0);
    vec3 Model_YN = vec3(0.0, -sign(Model_RayOrigin.y), 0.0);
    vec3 Model_ZN = vec3(0.0, 0.0, -sign(Model_RayOrigin.z));

As you can see, we are taking normals for the faces that go in the opposite direction of the ray. Those are the front faces.

Now let us figure out where the ray intersect the models. For that I'll be using the algebraic method for Ray-Plane Intersection, which looks like this:

    Ray   -> P = RayOrigin + t * RayDirection
    Plane -> P · N + d = 0
=>
    t = -(RayOrigin · N + d) / (V · N)
    P = RayOrigin + t * RayDirection

We, of course, do that for each of the three front planes, which looks like this:

    float Xd = -0.5;
    float Yd = -0.5;
    float Zd = -0.5;

    float Xt = -(dot(Model_RayOrigin, Model_XN) - Xd) / dot(Model_RayDirection, Model_XN);
    vec3 Model_PX = Model_RayOrigin + Xt * Model_RayDirection;

    float Yt = -(dot(Model_RayOrigin, Model_YN) - Yd) / dot(Model_RayDirection, Model_YN);
    vec3 Model_PY = Model_RayOrigin + Yt * Model_RayDirection;

    float Zt = -(dot(Model_RayOrigin, Model_ZN) - Zd) / dot(Model_RayDirection, Model_ZN);
    vec3 Model_PZ = Model_RayOrigin + Zt * Model_RayDirection;

---

This is where we get to the guesswork. I needed a way to check which plane is the correct one to use. My first guess was it would be the one further from the camera. And that is correct when we have all three front planes.

We can reuse Xt, Yt and Zt (which are distance along the ray) to figure out which point if further away from the camera. But we also need to discard the planes that we don't need. Since we want the further away, we can discard them by making the value 0. By trial and error I got it working like this:

    float Check_X = Xt * sign(floor(abs(Model_RayOrigin.x) * 2.0));
    float Check_Y = Yt * sign(floor(abs(Model_RayOrigin.y) * 2.0));
    float Check_Z = Zt * sign(floor(abs(Model_RayOrigin.z) * 2.0));

We, of course, still need to pick the one with greater value. Which will mean branching. If you can optimize the guesswork and the branching, kudos to you.

I also introduced a check to see if all these are greater than zero. If none of them are, it means the camera is inside the voxel volume.

    vec3 best = Model_BackFacePosition;
    if (Check_X > 0.0 || Check_Y > 0.0 || Check_Z > 0.0)
    {
        best = Model_PX;
        float best_length = Check_X;
        if (Check_Y > best_length)
        {
            best = Model_PY;
            best_length = Check_Y;
        }
        if (Check_Z > best_length)
        {
            best = Model_PZ;
        }
    }

Finally, you get your output:

    vec3 center_offset = vec3(0.5, 0.5, 0.5);
    vec3 ScaledPosition = (best + center_offset) * voxel_volume_size;
    vec4 o_Target = vec4(floor(ScaledPosition) / voxel_volume_size, 1.0);

---

At the end, I only used convert once. And the matrix I used with it was CameraToModel, which is inverse(Model) * inverse(View) (consider passing that as uniform). Replacing that into the code, this is the final code listing:

    mat4 CameraToModel = inverse(Model) * View;

    vec3 Camera_RayOrigin = vec3(0.0, 0.0, 0.0);
    vec3 Model_BackFacePosition = v_Position;
    vec3 Model_RayOrigin = (CameraToModel * vec4(Camera_RayOrigin, 1.0)).xyz;
    vec3 Model_RayDirection = normalize(Model_BackFacePosition - Model_RayOrigin);

    vec3 Model_XN = vec3(-sign(Model_RayOrigin.x), 0.0, 0.0);
    vec3 Model_YN = vec3(0.0, -sign(Model_RayOrigin.y), 0.0);
    vec3 Model_ZN = vec3(0.0, 0.0, -sign(Model_RayOrigin.z));

    float Xd = -0.5;
    float Yd = -0.5;
    float Zd = -0.5;

    float Xt = -(dot(Model_RayOrigin, Model_XN) - Xd) / dot(Model_RayDirection, Model_XN);
    vec3 Model_PX = Model_RayOrigin + Xt * Model_RayDirection;

    float Yt = -(dot(Model_RayOrigin, Model_YN) - Yd) / dot(Model_RayDirection, Model_YN);
    vec3 Model_PY = Model_RayOrigin + Yt * Model_RayDirection;

    float Zt = -(dot(Model_RayOrigin, Model_ZN) - Zd) / dot(Model_RayDirection, Model_ZN);
    vec3 Model_PZ = Model_RayOrigin + Zt * Model_RayDirection;

    float Check_X = Xt * sign(floor(abs(Model_RayOrigin.x) * 2.0));
    float Check_Y = Yt * sign(floor(abs(Model_RayOrigin.y) * 2.0));
    float Check_Z = Zt * sign(floor(abs(Model_RayOrigin.z) * 2.0));

    vec3 best = Model_BackFacePosition;
    if (Check_X > 0.0 || Check_Y > 0.0 || Check_Z > 0.0)
    {
        best = Model_PX;
        float best_length = Check_X;
        if (Check_Y > best_length)
        {
            best = Model_PY;
            best_length = Check_Y;
        }
        if (Check_Z > best_length)
        {
            best = Model_PZ;
        }
    }

    vec3 center_offset = vec3(0.5, 0.5, 0.5);
    vec3 ScaledPosition = (best + center_offset) * voxel_volume_size;
    vec4 o_Target = vec4(floor(ScaledPosition) / voxel_volume_size, 1.0);

---

Addendum: Optimizing

When we do dot product, we are taking the sum of the component wise product. Now look at how we defined Model_XN, Model_YN and Model_ZN. They only have one non-zero component. Thus we can replace the dot products with simple products on the corresponding components. While we are at it, let us make it a single vector Model_N.

We also have these Xd, Yd, and Zd. Make them into a vector d. Which also happens to be -center_offset.

Now, this code (I'm skipping Model_PX, Model_PY, and Model_PZ, soon we will not need them):

    vec3 Model_XN = vec3(-sign(Model_RayOrigin.x), 0.0, 0.0);
    vec3 Model_YN = vec3(0.0, -sign(Model_RayOrigin.y), 0.0);
    vec3 Model_ZN = vec3(0.0, 0.0, -sign(Model_RayOrigin.z));

    float Xd = -0.5;
    float Yd = -0.5;
    float Zd = -0.5;

    float Xt = -(dot(Model_RayOrigin, Model_XN) - Xd) / dot(Model_RayDirection, Model_XN);
    float Yt = -(dot(Model_RayOrigin, Model_YN) - Yd) / dot(Model_RayDirection, Model_YN);
    float Zt = -(dot(Model_RayOrigin, Model_ZN) - Zd) / dot(Model_RayDirection, Model_ZN);

Looks like this:

    vec3 Model_N = -sign(Model_RayOrigin);
    vec3 d = -center_offset;

    float Xt = -(Model_RayOrigin.x * Model_N.x - d.x) / (Model_RayDirection.x *  Model_N.x);
    float Yt = -(Model_RayOrigin.y * Model_N.y - d.y) / (Model_RayDirection.y *  Model_N.y);
    float Zt = -(Model_RayOrigin.z * Model_N.z - d.z) / (Model_RayDirection.z *  Model_N.z);

Well, make Xt, Yt and Zt a vector t, and we can write that like this:

    vec3 Model_N = -sign(Model_RayOrigin);
    vec3 d = -center_offset;
    vec3 t = -(Model_RayOrigin * Model_N - d) / (Model_RayDirection * Model_N);

We were using the old Xt, Yt and Zt to do so. With some added factors. Let us make those factors a vector too. So write this:

    vec3 f = sign(floor(abs(Model_RayOrigin) * 2.0));

Ok, those could be a vector too.. But wait! Why don't we just pick the best right away?

    float best_t = max(max(t.x * f.x, t.y * f.y), t.z * f.z);

Now the code inside the conditional can look like this (and we don't need Model_PX, Model_PY and Model_PZ anymore):

    best = Model_RayOrigin + best_t * Model_RayDirection;

This is how the code looks like now:

    mat4 CameraToModel = inverse(Model) * inverse(View);

    vec3 Camera_RayOrigin = vec3(0.0, 0.0, 0.0);
    vec3 Model_BackFacePosition = v_Position;
    vec3 Model_RayOrigin = (CameraToModel * vec4(Camera_RayOrigin, 1.0)).xyz;
    vec3 Model_RayDirection = normalize(Model_BackFacePosition - Model_RayOrigin);
    
    vec3 center_offset = vec3(0.5, 0.5, 0.5);

    vec3 Model_N = -sign(Model_RayOrigin);
    vec3 d = -center_offset;
    vec3 t = -(Model_RayOrigin * Model_N - d) / (Model_RayDirection * Model_N);
    vec3 f = sign(floor(abs(Model_RayOrigin) * 2.0));
    float best_t = max(max(t.x * f.x, t.y * f.y), t.z * f.z);
    vec3 best = Model_BackFacePosition;
    if (f.x > 0.0 || f.y > 0.0 || f.z > 0.0)
    {
        best = Model_RayOrigin + best_t * Model_RayDirection;
    }

    vec3 ScaledPosition = (best + center_offset) * voxel_volume_size;
    vec4 o_Target = vec4(floor(ScaledPosition) / voxel_volume_size, 1.0);

---

Godot

By the way, I used Godot to write the shader, so I could go straight to it (just add a mesh instance with a shader material, start typing and see the changes in real time). Godot shader language is based on GLSL, with some minor differences. This is the glue code I used:

shader_type spatial;
render_mode unshaded, cull_front;
uniform vec3 voxel_volume_size;

void fragment()
{
    mat4 View = INV_CAMERA_MATRIX;
    mat4 Model = WORLD_MATRIX;
    vec3 v_Position = (inverse(WORLD_MATRIX) * CAMERA_MATRIX * vec4(VERTEX.xyz, 1.0)).xyz;

    // THE CODE ABOVE HERE
    
    ALBEDO = o_Target.xyz;
}

And this is how it looks in engine:

Cheers!