Hyperbolic Voxel Sphere


I was trying to upload the voxel world to the gpu, and it worked (kinda) but meanwhile I found that the voxel world was not rendering correctly. Somewhere in my raytracer the rays are being bent, causing hyperbolic rendering instead of normal spatial rendering. This is extremely evident in the videos listed below:



I question why and how this bug exists but I will only understand once it is fixed.
Thanks for listening.

Comments

Log in with itch.io to leave a comment.

(8 edits)

Stating perhaps the obvious:  It looks to me like the ray origins and directions are not being calculated correctly.

Here are code snippets from my voxel ray march compute shader that calculate ray origins and directions (format is a little "wonky" pasted here in a code block) :

vec3 invViewProjOrigin(mat4 matrix) {
  vec4 origin = matrix * vec4(0, 0, -1, 0);
  return origin.xyz / origin.w;
}
vec2 ndcFragCoord(ivec2 resolution, ivec2 fragCoord) {
  return vec2((((fragCoord.x + 0.5) / resolution.x) * 2) - 1, 1 - (((fragCoord.y + 0.5) / resolution.y) * 2));
}
vec3 fragCoordRayDirection(ivec2 fragCoord, ivec2 resolution, float near, float far, mat4 invViewProjMatrix) {
  return normalize((invViewProjMatrix * vec4(ndcFragCoord(resolution, fragCoord), 1 + (near / far), 1)).xyz);
}
void main() {
  uint index = gl_GlobalInvocationID.x;
  if (index >= (resolution.x * resolution.y))
    return;
  ivec2 fragCoord = ivec2(index % resolution.x, index / resolution.x);
  vec3  origin    = invViewProjOrigin(invViewProjMatrix),
        direction = fragCoordRayDirection(fragCoord, resolution, settings.near, settings.far, invViewProjMatrix);
  ...
}

P.S.  I have my Vulkan pipeline's coordinates system y-flipped so positive y is up in world space.  Hence the "1 - (((fragCoord.y + 0.5) / resolution.y) * 2)" in ndcFragCoord().

Firstly, positive y should be up, anything else is heresy, and what is a projection matrix? I see everybody using one, but I don't use one myself. How do they work?

Google's AI describes it better than I can: "A projection matrix is a square matrix in linear algebra that, when multiplied with a vector, outputs the projection of that vector onto a specific subspace; essentially, it "flattens" the vector onto a lower-dimensional space while maintaining its relative position within that subspace."

I use perspective projection matrices for all things 3D in my voxel engine (calculations of ray origins and directions, rasterized lines and triangles).  Rasterizing pipelines in 3D graphics use the matrices to flatten projected 3D coordinates onto 2D screen space coordinates in order to rasterize in 2D.

(3 edits)

It is also strange that you say you think it is the ray origins and diretions that are not calculated correctly, (although you could very possibly be right, I wouldn't doubt it) But it is strange, because if I do a basic ray box intersection using the same ray directions and origins there is no warping at all. The warping exists entirely within the DDA algorithm I have. Here it is:


fn Sphere_DDA_March(
    start_pos: vec3<f32>,  // Starting position of the ray in world space
    end_pos: vec3<f32>,    // Ending position of the ray in world space
    voxel_size: vec3<f32>,       // Size of each voxel in the voxel grid
    bounds: vec3<i32>,     // Dimensions of the voxel grid (width, height, depth)
) -> vec4<f32> {
    var direction: vec3<f32> = normalize(vec3<f32>(end_pos - start_pos));
    var step: vec3<i32> = vec3<i32>(
    select(-1, 1, direction.x > 0.0),
    select(-1, 1, direction.y > 0.0),
    select(-1, 1, direction.z > 0.0)
    );
    var tMax: vec3<f32> = vec3<f32>(
        ((select(floor(start_pos.x / voxel_size.x) + 1.0, floor(start_pos.x / voxel_size.x), direction.x < 0.0) * voxel_size.x) - start_pos.x) / direction.x,
        ((select(floor(start_pos.y / voxel_size.y) + 1.0, floor(start_pos.y / voxel_size.y), direction.y < 0.0) * voxel_size.y) - start_pos.y) / direction.y,
        ((select(floor(start_pos.z / voxel_size.z) + 1.0, floor(start_pos.z / voxel_size.z), direction.z < 0.0) * voxel_size.z) - start_pos.z) / direction.z
    );
    var tDelta: vec3<f32> = vec3<f32>(
        abs(voxel_size.x / direction.x),
        abs(voxel_size.y / direction.y),
        abs(voxel_size.z / direction.z)
    );
    let center: vec3<f32> = vec3<f32>(15.0);
    var ray_world_coord = vec3<f32>(vec3<f32>(current_voxel) * voxel_size);
    var ray_dist: f32 = f32(sqrt(((ray_world_coord.x - camera.position.x) * (ray_world_coord.x - camera.position.x)) + ((ray_world_coord.y - camera.position.y) * (ray_world_coord.y - camera.position.y)) + ((ray_world_coord.z - camera.position.z) * (ray_world_coord.z - camera.position.z))));
    while(ray_dist < 256.0){
        ray_world_coord = vec3<f32>(vec3<f32>(current_voxel) * voxel_size);
        var world_cord = vec3<f32>(vec3<f32>(current_voxel) * voxel_size);
        ray_dist = f32(sqrt(((ray_world_coord.x - camera.position.x) * (ray_world_coord.x - camera.position.x)) + ((ray_world_coord.y - camera.position.y) * (ray_world_coord.y - camera.position.y)) + ((ray_world_coord.z - camera.position.z) * (ray_world_coord.z - camera.position.z))));
        var dist: f32 = f32(sqrt(((world_cord.x - center.x) * (world_cord.x - center.x)) + ((world_cord.y - center.y) * (world_cord.y - center.y)) + ((world_cord.z - center.z) * (world_cord.z - center.z))));
        if (dist < 60.0){
            let start_light: vec3<f32> = vec3<f32>(world_cord.x, world_cord.y + -0.5, world_cord.z);
            let end_light: vec3<f32> = vec3<f32>(world_cord.x, world_cord.y + 10.0, world_cord.z);
            let output = DDA_Light(start_light, end_light, voxel_size, bounds);
            if (output == 0) {
                return vec4<f32>(1.0, ray_dist / 100.0, 0.0, 1.0);
            } else {
                return vec4<f32>(0.75, ray_dist / 100.0, 0.0, 1.0);
            }
            
        }
        if (tMax.x < tMax.y && tMax.x < tMax.z) {
            current_voxel.x += step.x;
            tMax.x += tDelta.x;
        } else if (tMax.y < tMax.z) {
            current_voxel.y += step.y;
            tMax.y += tDelta.y;
        } else {
            current_voxel.z += step.z;
            tMax.z += tDelta.z;
        }
    }
    return vec4<f32>(0.0);
}

The code is too complicated for me to easily see where the calculations are incorrect.  Perhaps a test case in 2D would be easier to (visually) debug.

P.S. the bug is most prevalent at the poles of the sphere.
(this sphere has 6 poles where each face on a cube would be)